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Abstract. A summary of recent researches on nuclear dynamics with realistic microscopic quantum ap- 
proaches is presented. The Balian-Veneroni variational principle is used to derive the time-dependent 
Hartree-Fock (TDHF) equation describing the dynamics at the mean-field level, as well as an extension 
including small-amplitude quantum fluctuations which is equivalent to the time-dependent random-phase 
approximation (TDRPA) . Such formalisms as well as their practical implementation in the nuclear physics 
framework with modern three-dimensional codes are discussed. Recent applications to nuclear dynamics, 
from collective vibrations to heavy-ion collisions are presented. A particular attention is devoted to the 
interplay between collective motions and internal degrees of freedom. For instance, the harmonic nature of 
collective vibrations is questioned. Nuclei are also known to exhibit superfluidity due to pairing residual 
interaction. Extensions of the theoretical approach to study such pairing vibrations are now available. 
Large amplitude collective motions are investigated in the framework of heavy-ion collisions leading, for 
instance, to the formation of a compound system. How fusion is affected by the internal structure of the 
collision partners, such as their deformation, is discussed. Other mechanisms in competition with fusion, 
and responsible for the formation of fragments which differ from the entrance channel (transfer reactions, 
deep-inelastic collisions, and quasi-fission) are investigated. Finally, studies of actinide collisions forming, 
during very short times of few zeptoseconds, the heaviest nuclear systems available on Earth, are presented. 

PACS. 24.10.Cn Many-body theory - 21.60.Jz Nuclear Density Functional Theory and extensions (in- 
cludes Hartree-Fock and random-phase approximations) - 25.70.-z Low and intermediate energy heavy-ion 
reactions - 24.30.Cz Giant resonances 
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1 Introduction 

Nuclei are ideal to investigate fundamental aspects of the 
quantum many-body problem. They exhibit collective mo- 
tions built from coherent superpositions of the states of 
their constituents. Examples range from collective vibra- 
tions to the formation of a compound system in collisions. 
These features are common to other composite systems 
(atomic clusters, molecules, Bose-Einstein condensates...). 
Their study in nuclear systems is obviously part of a wider 
physics field. 

An interesting feature of the dynamics of quantum 
many-body systems is the interplay between collective mo- 
tion and single-particle degrees of freedom. The latter is 
a source of complex and fascinating behaviours. For in- 
stance, giant resonances are characterised by a collective 
vibration of many particles, but their decay may occur by 
the emission of a single nuclcon. Another example could 
be taken from the collision of composite systems where 
the transfer of few particles may have a strong impact on 
the final outcome of the reaction. 

To describe these complex systems, one needs to solve 
the quantum many-body problem. The description of the 
dynamics of composite systems can be very challenging, 
especially when two such systems interact. An important 
goal of nuclear physics is then to find a unified way to 
describe the dynamics of nuclear systems. Ultimately, the 
same theoretical model should be able to describe vibra- 
tion, rotation, fission, all the possible outcomes of heavy- 
ion collisions (elastic and inelastic scattering, particle trans- 
fer, fusion, and multifragmcntation), and even the dynam- 
ics of neutron star crusts. 

This desire for a global approach to nuclear dynam- 
ics has strongly influenced past research activities. Beside 
the quest for a unified model of nuclear dynamics, possi- 
ble applications of heavy-ion collisions to the formation of 
nuclear systems in extreme conditions (new exotic nuclei, 
large deformation and angular momentum, super-heavy 
elements...) are also strong motivations for the study of 
reaction mechanisms to optimise their production cross- 
sections. 

The Balian-Veneroni variational principle is a starting 
point leading to well-know formalisms, such as the time- 
dependent Hartree-Fock (TDHF) theory and the time- 
dependent random-phase approximation (TDRPA). Thanks 
to the recent increase of computational power, these ap- 
proaches have been applied to investigate several aspects 
of the nuclear dynamics with modern numerical codes. 
Recent applications to collective vibrations and heavy-ion 
collisions, as well as the underlying formalism and numer- 
ical details are discussed in this review article. 

The formalism and numerical developments are pre- 
sented in section 2. Studies of collective vibrations are 



2 The time-dependent Hartree-Fock theory 
and its extensions 

2.1 Introduction 

The quantum many-body problem is common to many 
theoretical fields [1,2]. It aims at describing the structure 
and dynamics of interacting particles. Electrons, atoms, 
molecules and nucleons are usual constituents of quantum 
many-body systems. 

In the non-relativistic regime, these systems obey the 
Schrodinger equation. 



' dt 



H\#(t)), 



(2.1) 



where |^(i)) describes the state of the many-body system 
and H is the microscopic Hamiltonian. We use the nota- 
tion h = 1. This equation can be solved exactly for simple 
cases only. Generally, one has to rely on some approxima- 
tions. 

Variational principles offer an elegant starting point to 
build such approximations. Indeed, they ensure an opti- 
mization of the equations of motion under the approx- 
imation that the variational space is limited to a sub- 
space of the full Hilbert (or Fock) space. Of course, with- 
out any restriction of the variational space, it is required 
that the chosen variational principle allows to recover the 
Schrodinger equation. However, their usefulness appears 
when restricting the variational space. Then, the validity 
of the approximation relies entirely on the choice of the 
variational space. In one hand, the latter has to be small 
enough so that the problem is numerically tractable. In 
the other hand, the variational space should contain the 
relevant degrees of freedom to allow for a realistic descrip- 
tion of the physical processes. 

Although the Schrodinger equation is unique, different 
variational principles have been developed in the past. One 
of the mostly used action is 



S[t ,ti;#(t)] = 



11 dt m)\(ij t -H 



The variational principle SS — is applied with the bound- 
ary conditions 5\P{t ) = S&ih) = 0. If 9 is allowed to 
span the entire Hilbert space, one recovers the Schrodinger 
equation. 

In most practical applications, mean-field models are 
considered in a first approximation, and, eventually, serve 
as a basis for beyond-mean-field approaches [3,4]. To con- 
struct such a mean-field theory from the above variational 
principle, one restricts the variational space by assuming 
that the N particles (we consider fermions) are indepen- 
dent. In this case, they may be described by a Slater de- 
terminant 



JV 



(2.3) 
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where dj creates a particle in the state \<fi) when it acts 
on the particle vacuum |— ). In such a state, all the in- 
formation is contained in the one-body density-matrix p 
associated to the single-particle operator 



JV 



(2.4) 



Solving the variational principle where the action defined 
in Eq. (2.2) is required to be stationary in the subspace 
of Slater determinants \<f>) with fixed boundary conditions 
at times to and t\ leads to (see appendix A) 



(2.5) 



where h[p] is the Hartree-Fock (HF) single-particle Hamil- 
tonian with matrix elements 



la/3 



and 



p a p = (<Pa\p\<fp) 



a' a a \ 



(2.6) 



(2.7) 



Eq. (2.5) is the time-dependent Hartree-Fock (TDHF) equa- 
tion. It was obtained by Dirac in 1930 [5]. It provides a 
self-consistent mean-field evolution where the interaction 
between the particles is replaced by a one-body mean- 
field potential generated by all the particles. It is, then, 
assumed that each particle evolves independently in this 
potential. 

We may question the validity of such a mean-field ap- 
proximation in the nuclear physics context. This assump- 
tion could only give an approximation of the exact dynam- 
ics and we do expect, in general, that the system will devi- 
ate from the independent particle picture. Indeed, the ex- 
act dynamics is given by the time-dependent Schrodinger 
equation [Eq. (2.1)], and, unless the Hamiltonian contains 
one-body operators only, the mean-field theory can only 
approximate the exact evolution of the system. The exact 
and mean-field Hamiltonians differ by the residual inter- 
action. The validity of the mean-field approximation de- 
pends, then, on the intensity of the residual interaction. 
The latter is a function of the state |<?) and, therefore, 
depends significantly on the physical situation. 

Starting from simple arguments [6], the time r over 
which the Slater determinant picture breaks down could 
be expressed as: 



~ 2 \N 



-1/2 



(2.8) 



a/3a/3 



where N is the number of particles, and a and (3 denote 
hole states while a and /3 denote particle states. In the 
nuclear physics context, typical values of the residual in- 
teraction lead to r ~ 100 — 200 fm/c. This gives an es- 
timate of the time during which one can safely consider 
that the independent particle approximation is valid. For 
longer times, as encountered, e.g., in heavy-ion collisions, 



the mean-field approximation can only be verified by com- 
parison with experiment and/or beyond mean- field calcu- 
lations. 

One important aspect of the TDHF approach is the 
treatment of one-body dissipation mechanisms. In partic- 
ular, one-body dissipation is crucial to properly describe 
low-energy heavy-ion collisions. Two kind of one-body dis- 
sipation mechanisms can easily be identified: 

— Coupling of collective motions with one-particle one- 
hole (Iplh) excitations. In the case of giant resonances, 
this coupling leads to the so-called Landau damping 
(see Sec. 3). In heavy- ion collisions, the field of the 
collision partner generates a distortion of the single- 
particle wave-functions. This is particularly true when 
part of the wave functions arc transferred from one 
fragment to another, leading to an excitation of the 
fragments, and, then, to a dissipation of the transla- 
tional kinetic energy. Another example is the case of 
fusion reactions. Multiple reflexions of single-particle 
wave-functions on the wall of the mean-field dissipate 
original collective translational kinetic energy into single- 
particle excitations and collective vibrations of the com- 
pound system. 

— Single-particle wave-function emission to the contin- 
uum. Emission of nuclcons into the continuum is a 
natural cooling mechanism of excited nuclei. As an ex- 
ample, the direct decay by nucleon emission induces 
an escape width of giant resonances (see Sec. 3.3). 

Naturally, extensions of TDHF including pairing and 
in-medium two-body correlations should be considered, at 
least to determine the validity of the mean-field approx- 
imation in terms of energy and simulation time. Larger 
variational spaces could then be considered in order to 
incorporate some effects of the residual interaction which 
are missing at the TDHF level. This is the case, for in- 
stance, with the inclusion of pairing correlations by taking 
a variational space of quasi-particle vacua, leading to the 
time-dependent Hartree-Fock-Bogoliubov (TDHFB) for- 
malism [7] . Another possible approach to include beyond 
mean-field correlations is to consider different levels of 
truncation of the Bogoliubov-Born-Green-Kirwood-Yvon 
(BBGKY) hierarchy [8-10]. The two first equations of this 
hierarchy read 



±Tr 

2 ±1 



2 [Wl2)Pl2j 



*^^Pl2 = [*1 + t 2 + ^Vu,pu] 



^Tr 3 [Vis 



^23,Pl23j , 

(2.9) 

where pi, pn and P123 are the one-, two-, and three-body 
density matrices, respectively We see that p\ and pn are 
coupled. In fact, each equation of the BBGKY hierarchy 
couples p\...k to p\...k+\ and forms a closed system of cou- 
pled equations equivalent to the Schrodingier equation for 
a finite many-body system. 

The two-body density matrix can be expressed as p\2 = 
PiP2(l — P12) + C12 where P\i is the permutation operator 
between particles 1 and 2, and C12 is the two-body cor- 
relation matrix describing, e.g., pairing correlations and 
in-medium two-body collisions. Neglecting C\i in the first 
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equation of the BBGKY hierarchy leads to the TDHF 
equation [Eq. (2.5)]. One could also include these two- 
body correlations by solving the coupled equations (2.9) 
and neglecting the three-body correlations. The resulting 
set of coupled equations for the evolution of pi and p\% 
is known as the time-dependent density-matrix (TDDM) 
formalism 1 . It is obvious that solving the TDDM equa- 
tions request much more computational efforts than stan- 
dard TDHF calculations. Only few applications have been 
made so far [11-17]. 

The difficulties of solving quantum many-body mod- 
els have sometimes been overcome by using their semi- 
classical limit. The semi-classical limit of the TDHF ap- 
proach is the Vlasov equation which is obtained by taking 
the Wigner transform of Eq. (2.5) and keeping only the 
first order in h. The in-medium two-body collision term 
can then be added, leading to the Landau- Vlasov equation 
which is the root of the Vlasov-Uehling-Uhlenbeck (VUU) 
[18] and Boltzmann-Uehling-Uhlenbeck (BUU) [19] trans- 
port theories. The latter approaches are common tools to 
interpret heavy-ion collision data at intermediate energy 
where the collision term is expected to play a significant 
role. However, at energies closer to the barrier, the colli- 
sion term is expected to be hindered by the Pauli block- 
ing. In addition, semi-classical approaches lead to a poor 
description of the ground state structure of the nuclei, 
such as their binding energies and deformations. As one 
of the goal of the present review article is to study the 
interplay between nuclear structure and dynamics, it is a 
natural choice to focus on fully quantum approaches only 
and semi-classical models won't be further discussed here. 

Several other variational principles have been proposed, 
depending on the physics one wants to address. For in- 
stance, Blaizot and Ripka introduced in 1981 a variational 
principle appropriate to the calculation of transition am- 
plitudes [7]. Inspired by this work, Balian and Veneroni 
(BV) presented the same year their variational principle 
for the expectation value of an observable [20] . In the lat- 
ter, both the state of the system and the observable of 
interest are allowed to vary in their own variational space. 

In particular, the BV variational principle allows a 
more transparent interpretation of the TDHF theory [20] : 
it is shown that the TDHF equation (2.5) is optimised for 
the expectation value of one-body observables. TDHF cal- 
culations are indeed successful to predict such quantities 
(several examples are presented in the following sections). 
It also explains why TDHF sometimes fail to reproduce 
other quantities such as fluctuations of one-body opera- 
tors [21-23], which are outside the variational space used 
to derive the TDHF equation. In fact, to predict such fluc- 
tuations, Balian and Veneroni proposed a prescription also 
based on their variational principle, but with a different 
variational space for the evolution of the observable [24, 
25]. 

In the following part of this section, we first describe 
the Balian- Veneroni variational principle. The latter is 
used to derive the TDHF equation as well as fluctua- 



1 Note that the TDDM formalism is not obtained by solving 
a variational principle. 



tions and correlations of one-body observables. Then, the 
Skyrme energy density functional describing the strong in- 
teraction between the nucleons is introduced. Numerical 
aspects are also discussed. Finally, perspectives for beyond 
TDHF calculations are presented. 



2.2 The Balian-Veneroni variational principle 

The BV variational principle has been applied to different 
problems in nuclear physics [26-32], hot Fermi gas [33], 
</> 4 -theory [34,35], and Boson systems [36,37]. The first 
realistic application of the BV prescription to fragment 
mass and charge distributions in heavy-ion collisions are 
reported in Ref. [32] and will be presented in section 4.4. 
The importance of the BV variational principle for the in- 
terpretation of the TDHF theory, which will be discussed 
and applied in the following sections, as well as the deriva- 
tion of the BV prescription for fluctuations and correla- 
tions of one-body observables, justify the more detailed 
discussion in this subsection. 

Let us first define two variational quantities: D(t), de- 
scribing the state of the system, and A(t), describing the 
evolution of the observable in the Hciscnbcrg picture. The 
application of the BV variational principle requires two 
boundary conditions: 



D(t ) = Do, 



(2.10) 



where the initial state of the system Dq is assumed to be 
known, and 



A{t l ) = A 1 , 



(2.11) 



where (Ai) is the final expectation value we want to com- 
pute at ti > to- 

The action-like quantity defined by Balian and Veneroni 
reads [20] 



J = Tr 

-/" 

J t 



A{ti)D{ti) 
dt Tr 



A(t)l^±+i{H(t),D(t)} 



. (2.12) 



We see that, imposing SaJ = 0, where 8a induces small 
variations of A(t), leads to the Liouvillc-Von Neumann 
equation 

.db(t) 



dt 



H,D(t) 



(2.13) 



which is fully equivalent to the Schrodinger equation. 

To get Eq. (2.13), we first note that, according to the 
boundary condition in Eq. (2.11), A(ti) is fixed and we 
get 



SaJ 



dt Tr 



to 



S A A(t) t^L +i [H(t),D(t)] 



To be equal to zero for any variation of A, the term inside 
the brackets must be zero, giving Eq. (2.13). 
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Variations of D(t) should also be considered. It is easier 
to first rewrite Eq. (2.12) as 



J = Tr 



A(t )D(t ) 



dt Tr 



D(t) (^-+i[H(t),A{t)) 



(2.14) 



To get Eq. (2.14) we integrate by part the / dtAdtD term 
in Eq. (2.12) and we use the relation 

Tr(A[H,b}) = -Tt(D[H,A}). 



Equivalently to Eq. (2.13), with the boundary condi- 
tion in Eq. (2.10), requiring 8dJ = leads to 



dA(t) 
1 dt 



H, Alt) 



(2.15) 



which is also equivalent to the Schrodingcr equation, and 
is expressed in the Heiscnberg picture. 



2.3 Derivation of the time-dependent Hartree-Fock 
equation 

The TDHF theory is obtained under the approximation 
that A(t) is constrained to be a one-body operator and 
that D(t) is an independent particle state for all t. As a re- 
sult, TDHF is optimised for the expectation value of one- 
body operators. On the contrary, other quantities, such as 
expectation values of two-body operators, may be not well 
predicted because they are outside the variational space 
of A{t). 

To get the TDHF equation, we impose that the vari- 
ation 8a leaves A in the space of one-body operators. As 
we consider arbitrary variations, we can choose 



SAA(t) = d^d/3 for t < t < t\ 



(2.16) 



and SA(ti) = due to the boundary contition in Eq. (2.11). 
Requiring Sa J = 0, we get from Eq. (2.12) 



Tr 



dD 
~~dt 



4«/3 I ^jr+i[H,D] 



= 0. 



(2.17) 



In addition, the state of the system is constrained to be 
an independent particle state. The variational space for 
D(t) is then defined by D(t) = \(f>{t)) (<f>(t)\ where \<f>(t)) is 
a Slater determinant. Using the one-body density matrix 
defined in Eq. (2.7), we get 



d\dfs,H 



Consider a Hamiltonian of the form 



-ijki a\d^aiak 

ijkl 



(2.18) 



(2.19) 



where matrix elements associated to the kinetic energy 
and to the anti-symmetric two-body interaction are given, 
respectively, by 



tij = — (i\f\j) and 

Vijkl = Vijkl — Vijik- 



(2.20) 
(2.21) 



Reporting the Hamiltonian expression [Eq. (2.19)] in 
Eq. (2.18), we get 



7— 

dt 



+ 4 X! Vki m n([alap,alal a n a m ]), (2.22) 

klrrtn 

where the time variable has been omitted and (• • ■) denotes 
the expectation value on \4>(t)) to simplify the notation. 
Eq. (2.22) leads to the TDHF equation 



(2.23) 



The single-particle Hartree-Fock Hamiltonian reads 

h[p]=t + U[p] (2.24) 
with the self-consistent mean-field 

U[p]ij = ^ikjl Plk- (2.25) 

kl 

To show the equivalence between Eq. (2.22) and the 
TDHF equation, let us start with the term associated to 
the kinetic energy: 

( [ 4 o-i . 4 k ] ) = s jk (4 k) - (4 4 o-j h i) 
-5u (4 flj) + (4 4 aidj) 
= Sjk pu — <\i Pjk- 

The kinetic energy term reduces to 

2J t k i([a\ aj , a[ai}) =Y^ (tjk Pkt - t ki p jk ) = [t, p]^ . 

kl k 

For the two-body interaction, we need the expectation 
value of the commutator 



[a\a,j , 4 4 a«a m " 



(4 4 a 



«n 0,m) 3jk ~ ( a i &fc On &m) $jl 



+ (4 4 4 dj d n d m ) - (4 4 ^" 



+ {d\ dj d m a,j) S nl - (aj. a[ d\ d n d m 

The two terms with 6 annihilation/creation operators can- 
cel out and we get 

([4«J . 44«««™]> 

= {PmiPnl — PmlPni) &jk + (PmkPni ~ PrniPnk) Sjl 
+ (PjlPnk — PjkPnl) &mi + {PjkPml ~ PjlPmk) 4ii 
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Altogether, the two-body interaction contribution reduces 
to 

4 Vkimni [a\ aj,&l aj a n a m ] ) 

klmn 

= ^ [VjMm (PtiPmk ~ PlkPmi) 
klm 

+Vklim {PjlPmk ~ PjkPml)] 

= ^ [Vjklm (piiPmk) — Vklim (PjkPml)] 
klm 

= ( U lP\jkPM - U[p]kipjk) = [U[p\iP] J% 

k 

where we have used v klmn = ~v k i nm = —Vikmn- Gathering 
the kinetic and interaction terms gives Eq. (2.23). 

We see that the TDHF equation is obtained by solving 
the BV variational principle with the variational spaces 
restricted to Slater determinants for the state of the sys- 
tem and to one-body operators for the observable. It is 
interesting to see that only the variation of A is needed to 
get the TDHF equation. In addition, solving the TDHF 
equation allows to compute any one-body observable, and 
the equation does not depend on the final time t\ entering 
the BV action. These properties are specific to the TDHF 
case. In a more general case, the resulting equations of 
motion are obtained from both the variation of A and D, 
and the results are valid for only one observable A\ and 
one final time t\. 

We recall that the TDHF equation is optimised for 
the expectation value of one-body operators and its pre- 
dictive power for other purposes may be questionable. In 
particular, two-body operators and fluctuations of one- 
body operators are outside the variational space of the 
observable used to derive the THDF equation. In the next 
section, we solve the BV variational principle in order to 
compute such fluctuations. (The case of more general two- 
body operators would be more complicated and the re- 
sulting equations of motions are not expected to be easily 
solvable numerically.) 



2.4 Fluctuations and correlations of one-body 
observables 

Let X and Y be two observables. Their correlation in the 
state 1^) is defined as 

a XY = ^((XY) + (YX))-(X)(Y). (2.26) 

The case X = Y defines the fluctuations of X: 

a xx = a x = ^{X^)-{X)\ (2.27) 

If X and Y are one-body operators, we see that o\ Y 
includes the expectation value of the square of a one-body 



operator, which contains a two-body contribution. 
Indeed, 

XY = 2J X a pa\flp ^ Y^a^a^ 

aft iiv 

The last term is clearly of a two-body nature. The TDHF 
theory is then not optimised for the prediction of correla- 
tions and fluctuations of one-body operators. 

A possible improvement would be to solve the BV 
variational principle with a variational space for the ob- 
servable which includes square of one-body operators [38] . 
This approach, however, leads to complicated equations of 
motion due to an intricate coupling between the evolution 
of the observable A{t) and of the state D(t). 

Fluctuations and correlations of one-body operators 
Qi = 'Yl, ol aQi a p<i\ l Q>p can also be computed from the ex- 
pectation value of 

Ax = e-^ ei ^ (2.28) 
in the limit Si — > 0. Indeed, 

ln(i0 = -]T><Qi) + \Y, £ i e 3dj +0(e 3 ) 1 (2.29) 

i " ij 

where Cy = Oq . q . . The linear and quadratic dependences 
in e of ln(.Ai) then lead to the expectation values and 
fluctuations/correlations of the one-body observables Qi, 
respectively. 

In Ref. [24], Balian and Veneroni applied their varia- 
tional principle with a variational space for the observable 
parametrised by 

A(t) = e- L{t) = e-Ea^Wata^ ( 2 .30) 

while the density matrix is constrained to be an indepen- 
dent particle state. For a pure state, the latter takes the 
form 

b(t) = \m(Ht)\, (2.3i) 

where \4>(t)) is a Slater determinant. In fact, the origi- 
nal derivation of Balian and Veneroni [25] involves more 
general mean-field states of the form 

D(t) = e -m(t)-M(t) _ e -m(t)-J2 afl Mcpitjalap ^ (2.32) 

The particular case of a Slater determinant can be ob- 
tained with this paramctrisation by letting the eigenvalues 
of the matrix M{t) tend to ±00, with the normalisation 
factor m(t) also tending to +00 in such a way that the 
norm z(t) = Ti(D{t)) is equal to 1 [39]. In this case, the 
Slater determinant \4>{t)) is built from the eigenvectors 
of M(t) associated with the eigenvalues —00. Compar- 
ing Eqs. (2.30) and (2.32), we see that both variational 
spaces for the state and the observable are similar, being 
both composed of exponentials of one-body operators. 
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The resolution of the BV variational principle with this 
choice of variational spaces and assuming the boundary 
conditions given in Eqs. (2.10) and (2.11) can be found in 
Refs. [25,31] and in Appendix B. 

Expanding the one-body density matrix in terms of e, 



p(t)=pW(t) + pW(t) + 0(e 2 ), 

where p^> is of the order e, the main results read: 
— The expectation value of Qi is given by 

(Qi)(ti)=tr(pW(t 1 )Q i 



(2.33) 



(2.34) 



where p(°> is given by the TDHF equation (2.5) with 
the boundary condition 



pi>o)=Tr 



(2.35) 



and Do is the initial density matrix. 
— The fluctuations/correlations OQ i Q i 



TDHF 



J 

where 




(% DHF (ti) 



ij 



P <- Hti)QjO--p ( -°Hh)) 



2 
tr 



are the (square of the) fluctuations/correlations ob- 
tained from the standard TDHF approach. 

Eq. (2.36) gives fluctuations and correlations which differ 
from the standard TDHF result [Eq. (2.37)]. This is not 
surprising as only Eq. (2.36) is optimised for these fluctu- 
ations/correlations. The result in Eq. (2.36) takes into ac- 
count possible fluctuations around the TDHF mean-field 
evolution in the small amplitude limit, i.e., at the RPA 
level [24,25,40]. 

The additional term in Eq. (2.36) involves p^\ i.e., 
the part of the one-body density matrix which is linear 
in s. However, the equation of motion for the latter is not 
trivial. Fortunately, for a Slater determinant, Eq. (2.36) 
can be re- written so that it is easier to implement. The 
final result reads (see appendix B) 



1 



lim 

ei.Ej-s-O 2Si£j 



( 



P {0) (t )-m(to,£ l )) 
p (0) (to) -»&(*«>,£,•))] , (2-38) 



where the single-particle matrices r](t,e) obey the TDHF 
equation (2.5) with a boundary condition defined at the 
final time t\. 



The fluctuations <jQ i 
Qi = Qj , leading to 



are determined by taking 



Cu{ti) 



lim — ?rtr 

e,^o 2e? 



(V >(to)-^(io,e*))' 



(2.40) 



(2.39) 



with the boundary condition given in Eq. (2.39). 

The fluctuations are generated by the boost in Eq. (2.39) 
and propagated in the backward Hciscnberg picture from 
t\ to t according to the dual of the time-dependent RPA 
equation. This is why Cy(ti) is expressed as a function of 
density matrices at the initial time to. It is easy to show 
that, if the backward trajectories rji have the same mean- 
field as the forward evolution, then Eq. (2.38) leads to the 
TDHF expression in Eq. (2.37). If, however, deviations oc- 
cur around the original mean-field, then additional terms 
appear and lead to an increase of Cij(t\). 

Eq. (2.39) imposes to solve the TDHF equation first 
for p(°\t) forward in time, and then for r](t) backward in 
time. Numerical applications solving Eq. (2.40) with the 
boundary condition given in Eq. (2.39) are detailed in [26- 
28,30-32]. In practice, several (typically ~ 5 in [32]) back- 
ward TDHF trajectories with different but small values of 
e are performed to compute the limit in Eq. (2.40). Fig. 2.1 
gives a schematic illustration of the numerical technique. 
In Ref. [32], Eq. (2.38) is also solved. It is used to compute 
the correlations between proton and neutron numbers in 
fragments following deep-inelastic collisions in addition to 
their fluctuations (see section 4.4). 

It is important to note that, although Eqs. (2.38) and 
(2.40) are rather easy to solve numerically with a slightly 
modified TDHF code, they only provide the fluctuations 
and correlations for a specific choice of Qi and Qj at the 
time t\. Another choice of operator(s) and/or final time 
implies to compute numerically another set of backward 
TDHF evolutions. This was not the case for the calculation 
of expectation values of one-body observables which only 
requires one forward TDHF evolution. Indeed, solving the 
TDHF equation is an initial value problem and the result- 
ing one-body density matrix can be used to compute any 
one-body observable (i.e., the TDHF equation depends 
neither on t\ nor on Qi). As a result, the determination 
of fluctuations and correlations of one-body observables 
from the BV variational principle is much more computa- 
tional time consuming than standard TDHF calculations 
of their expectation values. 

Finally, it is worth mentioning a recent alternative 
derivation of Eq. (2.38) considering small fluctuations in 
the initial state at time to and propagated within the 
stochastic mean-field (SMF) approach [40] . Numerical ap- 
plications of the SMF method can be found in Refs. [41- 
43]. 



2.5 Skyrme Energy-Density-Functional 

Early TDHF calculations [44,45] were based on simpli- 
fied Skyrme effective interactions [46]. Modern applica- 
tions, however, use full Skyrme functionals including spin- 
orbit interactions to reach a realistic description of the 
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where 6 = S(r(l) - f (2)), k = ^ (V(l) - V(2)) (relative 
momentum), k' complex conjugated of k acting on the left, 
and R = (f (1) + r(2)) /2. The operators & = &(1) + <r(2), 
with = & x (i) e x + a y (i) e y + a z (i) e z , are expressed in 
terms of the Pauli matrices o x / y / z {i) acting on the spin of 

the particle i. P a = [1 + <x(l) ■ ct(2)] /2 corresponds to the 
exchange of the spins. The particle density in r is noted 
P( r ) = P( rs 9i rs l) where p is the one-body density 
matrix, s the spin and q the isospin. The "ti" and "t 2 " 
terms are non-local in space and simulate the short range 
part of the interaction. Finally the last term accounts for 
the spin-orbit interaction. 

The very interesting aspect of this interaction is its 
zero range nature, which greatly simplifies the mean-field 
expression in coordinate space (see below). Parameters 
(to-3, £0-3! Wq and a) are generally adjusted to reproduce 
nuclear properties like saturation and incompressibility of 
nuclear matter and selected properties of finite nuclei (see 
for instance [49]). 

In the following of this section, when there is no ambi- 
guity, the spin and isospin indices, as well as the indices a 
denoting single-particle states are omitted in the notation 
for clarity. We also assume that p is diagonal in isospin. 
Let us define the usual densities entering the Skyrmc EDF: 



Fig. 2.1. Schematic illustration of the main numerical steps 
to compute the fluctuation of the one-body observable N (e.g., 
neutron numbers in one fragment following a heavy-ion colli- 
sion) . (a) A standard TDHF evolution is performed from to to 
ti. (b) The transformation given in Eq. (2.39) is applied on the 
one-body density matrix at time ti with several (small) values 
of s. (c) Backward TDHF evolutions from ti to to are per- 
formed for each e and the resulting one-body density matrices 
at time to are used to compute the (square of the) fluctuations 
from Eq. (2.40). 



strong interactions between the nucleons. Recent studies 
with the Gogny effective interaction [47] have also been 
reported [48]. 

Modern Skyrme energy density functionals (EDF) are 
usually derived from a Skyrmc effective interaction of the 
form 



6(1, 2)= to [l + xoPa) 5 

+ i h (i + xi p ff ) (k 2 s + s k 2 

+ i 2 (l + x 2 P a ) (k'-Sk 
1 



6 



t 3 1 + x 3 P a P a (R) 5 



+ iW or ■ fk' x <5k 



as as 

r(r)=J2 IV^I 2 (2.43) 

as 



Z 'i 

V.J(r) = V99*(rs) x V</?(rs') ■ (s\a\s')(2A5) 

ass' 

S(r)=]T ^(r S )^(r.s') (s\cr\s'), (2.46) 



where c.c. means "complex conjugated", are the local, ki- 
netic, current, (gradient of) spin-orbit, and spin densities, 
respectively. The j and S densities are time-odd and vanish 
in time-reversal invariant systems. They are, however, im- 
portant in time-dependent calculations to ensure Galilean 
invariance [50]. 

The total energy of an interacting system can be writ- 
ten as an integral of a local energy density 



(2.41) 



E = J dvH{r) 



(2.47) 
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Within the framework of the Skyrme EDF, we have [51] 
H(r) = £ T + Bl/?2 + i?2 ^ 

9 

+B 3 (pr - j 2 ) + Bi J2 (Pi T i - f q ) 

9 

+B 5 pAp + B 6 Y, PgM + B 7P 2+a + B sP a p\ 
i i 

+B 9 (pV-3 + j Vx S + ^ Pq V-3 q + i g - Vx S q ) 

q 

+B 10 S 2 + B u Y, s2 q + Bi2P a S 2 + B 13P a ]T S 2 , 

9 9 

(2.48) 

where the coefficients Bi are the usual Skyrme functional 
coefficients: 



B x 

B 2 

B a 
Bi 

B 5 
B 6 
B 7 
Bs 

B 9 
Bio 
B n 

Bl2 
B13 



to 



1 



2 



to 



x 



1 + i 



t2 1 + 4 



1 

"4 
1 

16 



h xi + 



*2 X 2 



1 + y) 



1 

To 

h 

12 



3*1 ( Xl + - ) + t 2 | X2 



1 

2 

X2 

Y 
i 

' 2 



t 3 



1, 



2 

tpxp 
4 

J_o 

4 

24 
24' 



Wo 



(2.49) 



In fact, the complete Skyrme functional is more gen- 
eral and contains other terms of the form S • AS and 
with other densities, i.e., the spin-current pseudo-tensor 

J and the spin kinetic energy density T [50,52]. All or 
some of these additional terms are sometimes included in 
TDHF calculations [52-54]. In Eq. (2.48), only the anti- 

symmetric part of J, which is the spin-orbit density J, is 
included. The spin-orbit energy is indeed expected to be 
more important (by about one order of magnitude) than 
the other spin-gradient terms [49]. 



The Skyrmc-HF mean-field is derived from Eq. (2.6) by 
replacing the expectation value of the Hamiltonian on the 
Slater determinant by the expression of the Skyrme EDF. 
The action of this field on single-particle wave functions 
is then given by [51] 

(h[p](p Q ) (r, s) = 
_ h 2 



E 



V+C/Jr)+iCJr)-VU, 



2S Q (r) 

+V g „(r) • {s\<T\s')+iW ga (r) ■ ((s\<r\s') x V) ] tpjr^'). 

(2.50) 

The derivatives act on each term sitting on their right, 
including the wave function. The fields (functions of r) 
read 



2m* q 



2m 



B 3 p + B± p q 



(2.51) 



U„ = 



2B lP + 2B 2 p q + B 3 (T + iV-j) + B±(T q + iV • j,) 
+2B 5 Ap + 2B 6 A Pq + (2 + a)B 7 p 1+a 

+B 8 [ap a - 1 Y,p 2 q + 2p a P q \ + B 9 (V ■ J + V ■ J q ) 

q 

+ap a - 1 (B 12 S 2 +B 13 Y,S 2 q ) (2-52) 

i 

B 9 V x (j+j (Z ) + 2B 10 S + 2B 11 S f; 
+2p a {B 12 S + B 13 S q ) (2.53) 
-B 9 V(p + Pq ) (2.54) 
2B 3 j + 2B 4 j,-B9Vx(S + S ? ), (2.55) 



where the derivatives act on the first term sitting on their 
right only. 

In addition to this mean-field potential, protons are 
also affected by the Coulomb interaction. The direct part 
of the Coulomb energy reads 



E. 



<lrr 



— I ctr 



d 3 r , Pp( r )P P ( r ') 



(2.56) 



The latter is usually computed by solving first the Poisson 
equation to get the Coulomb potential V c (r), and, then, 
by evaluating the integral i J d 3 rp p V c . The exchange part 
of the Coulomb energy is usually determined within the 
Slater approximation as 



Ei 



-3e 2 fi 

IT 



d 3 rp p (r)3. 



(2.57) 



As a result, the contribution of the Coulomb interaction 
to the proton mean-field reads 



U c = Vc 



3pp 
n 



(2.58) 



2.6 Numerical solution of the TDHF equation 

As we saw in section 2.4, the calculations of both expec- 
tation values and fluctuations/ correlations of one-body 
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obscrvablcs imply to determine the time evolution of the 
one-body density matrix with the TDHF equation (2.5). 
Few numerical codes solving the TDHF equation in three 
dimensions with a full Skyrme EDF including spin-orbit 
terms [46,49] are now available [55-59]. They were used 
in many studies of giant resonance properties [60-62,57, 
58,63,30,64-66], neutron star crust properties [59,67] and 
heavy-ion collisions. Amongst the latter, fusion and nucleus- 
nucleus potentials [55, 68, 69, 53, 52, 70-79, 41, 80-91] , isospin 
equilibration [68,74,75,92-95,89,32,96], Coulomb excita- 
tion [69,71,72,75], rotational properties [97], quasi-elastic 
transfer [76, 98, 42, 99, 100, 43, 96], breakup [17], deep-inelastic 
collisions [32], as well as the dynamics of a-clustcrs [101] 
and actinide collisions [102,103] have been investigated. 



of mass motion, i.e., by replacing the total kinetic energy 
f by 

P 2 A (Ef=l Pi^. 



T 



2Am 



T 



T 



2Ar 



2Ar 



p* p j 



(2.59) 



2.6.1 Numerical method 

Most of these applications were performed on a three di- 
mensional cartesian grid using a time iterative method. 
Wc give here a summary of the main steps which are usu- 
ally followed to treat the collision of two nuclei: 

1. Static Hartree-Fock (HF) calculations are performed 
to determine the initial condition where the nuclei are 
usually assumed to be in their HF ground state. 

2. The nuclei are placed in a larger box, avoiding any 
overlap of the HF solutions. This latter condition al- 
lows to construct a single Slater determinant from the 
two initial HF states. 

3. A Galilean boost 2 is applied at the initial time assum- 
ing that the nuclei followed a Rutherford trajectory 
prior to this time. 

4. The TDHF equation is solved iteratively in time and 
expectation values of one-body observables are even- 
tually computed at each time step to get their time 
evolution. 

Of course, variations of the main numerical steps described 
above are possible. For instance, one can question the va- 
lidity of the assumption that the nuclei are in their HF 
ground state at initial time. In particular, heavy nuclei 
generate strong Coulomb fields which may induce long 
range excitations of the collision partners [69] . 



2.6.2 Center of mass corrections in heavy-ion collisions 

HF calculations arc usually performed with center of mass 
corrections to improve the description of the nucleus in its 
intrinsic frame. This is done by removing spurious center 



2 In case of a single nucleus, for instance to study its response 
to a specific excitation, the Galilean boost is replaced by the 
appropriate velocity boost generating the excitation (exam- 
ples are given in section 3). Alternatively, one can start with 
a constrained Hartree-Fock (CHF) solution obtained with an 
external constraint in the HF calculation. The response to the 
excitation is then studied by relaxing the constraint in the 
TDHF calculation [104]. 



This correction contains a one-body and a two-body con- 
tributions. Usually, the two-body term is neglected and 
only the one-body part of the correction is included in 
standard HF calculations [49] . 

The initial condition of a TDHF calculation uses static 
Hartree-Fock (HF) or constrained Hartree-Fock (CHF) so- 
lutions. However, center of mass corrections are difficult 
to incorporate in TDHF calculations of heavy-ion colli- 
sions [81]. They are usually neglected to allow a consistent 
treatment of colliding partners. Indeed, these corrections 
are explicitly dependent on the number of nucleons Aj of 
the collision partner i and would induce a different treat- 
ment of the single-particle wave-functions depending on 
which nucleus they come from. 

To treat structure and dynamics on the same footing, 
one can also neglect the center-of mass corrections in the 
initial HF calculations [55]. However, one should then use 
an EDF which has been fitted without these corrections. 
This is the case of the SLy4d parametrisation [55] of the 
Skyrme EDF which is widely used in this review article. 



2.6.3 Numerical approximations and algorithm 

Most of the numerical approximations and techniques used 
in modern TDHF codes are based on similar algorithms, 
but may contain differences such as in the calculation of 
spatial derivatives. For instance, the latter are computed 
with finite-difference formulae in the tdhf3d code [55], 
while spline and fast-Fourier transform (FFT) techniques 
are used in Refs. [58] and [56], respectively. Typical regu- 
lar mesh spacing with Ax ~ 0.6 fm [64], 0.8 fm [55] and 
1.0 fm [58,56] are used, although adaptive grids have been 
also considered [57]. 

The tdhf3d code contains all the time-odd and even 
terms of standard Skyrme EDF (see Sec. 2.5). The inclu- 
sion of time-odd terms is indeed crucial for a proper de- 
scription of translational motion and to avoid spurious ex- 
citations [53]. As time-reversal symmetry is not assumed, 
the code contains no degeneracy of single-particle wave- 
functions. This means that up to ~ 500 (for, e.g., two 
actinides) wave- functions are evolved in time. 

Due to the self-consistency of the mean-field, the TDHF 
equation needs to be solved iteratively in time with a small 
time step increment At with typical values ranging from 
- 5xl0~ 25 s [57,64] to~ 1.5xl0~ 25 s [55,58]. Over small 
time intervals [t,t + At], the HF Hamiltonian is assumed 
to be constant. However, to conserve energy, the numeri- 
cal algorithm should be symmetric with respect to time- 
reversal operation. This implies to consider the Hamilto- 
nian at time t + 4f- for the evolution of single-particle 
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Fig. 2.2. (courtesy of D. Lacroix). Illustration of the different 
kinds of correlations. From left to right: in medium particle- 
particle collisions, pairing correlations, and large amplitude dy- 
namical correlations. 



wave-functions from t to t + At [44] 

\v(t + At)) we-'TK'+f) |„(t)). 



(2.60) 



A schematic illustration of the real time propagation could 
be written as: 



(") 



(n) 



>} 



P 



(n) 



fr 



,(»+*) 



|^+ 1 )) = e -^^V| n) 



(2.61) 

where |z/ n )) corresponds to an approximation of \u(t n = 
nAt)). In this algorithm, starting from the density at time 
t, a first estimate of the density at time t + At, denoted by 
p(n+i) j g kt amec L The Hamiltonian used in the propaga- 
tor is computed using the average density obtained from 
and p( n+1 \ Then, the new density at time t + At is 
obtained using this Hamiltonian. An approximate form of 
the exponential function is generally used which in some 
cases, breaks the unitarity, and orthonormalisation of the 
single particle states must be controlled. 



2.7 Perspectives for beyond TDHF calculations 

The TDHF theory is a microscopic quantum mean-field 
approach of independent particles. Extensions of the TDHF 
theory include correlations which are not present at the 
TDHF level. These correlations are typically of three kinds 
illustrated in Fig. 2.2: 

— In medium particle-particle collisions. In a first ap- 
proximation, this collision term can be neglected thanks 
to the Pauli blocking. However, in violent collisions, 
where the Pauli principle is less efficient to block colli- 
sions between nucleons of the two colliding fragments, 
this term is expected to affect the dynamics. It is also 
responsible for the thermalisation of the compound nu- 
cleus and for the spreading width of giant resonances. 

— Pairing correlations. They generate a superfluid phase 
in nuclei. They have a strong effect on the structure of 
mid-shell nuclei (e.g., odd-even mass staggering) and 
in transfer reactions where they favour the transfer of 
paired nucleons. 



— Large amplitude dynamical correlations. Unlike a sin- 
gle Slater determinant which is usually localised in a 
potential energy surface (PES), the correlated state 
may be described by a configuration mixing of localised 
states across the entire PES. A typical example is the 
zero point motion along a collective coordinate. These 
correlations allow also for the state to be in a classically 
forbidden region of the PES, and, then, are necessary 
to treat quantum tunnelling of the many-body wave 
function, e.g. in sub-barrier fusion. 

The Balian-Veneroni prescription discussed before is 
an example where fluctuations and correlations of specific 
observables, not included at the TDHF level, arc described 
dynamically. In the classification used above, the included 
correlations belong to the class of dynamical correlations, 
although, in this case, they are of small amplitude na- 
ture. In fact, the fluctuations obtained in Eq (2.40) arc 
those included in the time-dependent random phase ap- 
proximation (TDRPA) which is obtained assuming small 
fluctuations of the density matrix around the average evo- 
lution [25,40]. 

A possible extension of TDHF including dynamical 
correlations of large amplitude involves the path integral 
technique with the stationary phase approximation (SPA) 
[45] . The SPA assumes that the path integral is dominated 
by the classical action (for a single particle), or, equiva- 
lently, by the TDHF action for a many-fermion system. 
The method provides an elegant way to include fluctu- 
ations around the mean-field trajectory. Unlike the BV 
prescription, however, these fluctuations are not limited 
to the small amplitude limit. This approach leads to a 
self-consistent eigenvalue problem in four space-time di- 
mensions. It is in fact similar to a HF eigenvalue problem 
with time as an additional dimension. In case of a vibra- 
tional motion, the problem involves the (difficult) task to 
find periodic solutions of the TDHF equations. Another 
important possible application of this path integral ap- 
proach is to treat quantum tunnelling of the many-body 
wave function through a barrier. In this case, quantum 
many-particle closed trajectories in imaginary time need 
to be found. Up to now, realistic applications in imagi- 
nary time have faced the difficulties brought by the lim- 
itations of computational power, and only simple cases, 
such as the spontaneous fission of the 8 Be in two a have 
been studied [45]. However, the recent increase of compu- 
tational power, and the strong interests in understanding 
dcep-sub-barrier fusion [105], should lead to a revival of 
these techniques. 

An alternative way to include large amplitude dynam- 
ical correlations is to derive collective Hamiltonians, such 
as the Bohr Hamiltonian, from microscopic calculations 
using the generator-coordinate method (GCM) [106,107] 
and its time-dependent extension (TDGCM) [108]. Recent 
applications of the GCM to low-lying collective excita- 
tions can be found in Refs. [109-111]. The TDGCM has 
also been applied with the Gaussian overlap approxima- 
tion (GOA) to investigate the fission process [112]. 

Other examples of beyond-TDHF approaches include 
the time-dependent Hartree-Fock-Bogoliubov (TDHFB) 



12 



Cedric Simenel: Nuclear Quantum Many-Body Dynamics 



theory [2,113,114,65], the time-dependent density-matrix 
(TDDM) theory [11-17], the extended TDHF [115-120] 
and stochastic TDHF theories [121-123,4], and stochas- 
tic mean-field based on functional integrals [124,125]. Ta- 
ble 2.1 gives a summary of these approaches. Numerical 
applications of the TDHF, BV, and TDHFB approaches 
are presented in this review. 



3 Collective vibrations 

3.1 Introduction 

A particular interest in strongly interacting systems is 
their ability to present disorder or chaos, and, in the same 
excitation energy range, well-organised motion. Atomic 
nuclei are known to show both behaviors [126]. In partic- 
ular, they exhibit a large variety of vibrations, from low- 
lying collective modes to giant resonances (GR) with exci- 
tation energy usually above the particle emission thresh- 
old [127]. 

Baldwin and Klaiber observed the isovector giant dipolc 
resonance (GDR) in photofission of uranium nuclei [128], 
interpreted as a vibration of neutrons against protons [129] . 
Other kinds of GR have been discovered, such as the 
isoscalar giant quadrupole resonance (GQR) associated 
with an oscillation of the shape between a prolate and an 
oblate deformation [130], and the isoscalar giant monopolc 
resonance (GMR) corresponding to a breathing mode [131— 
133]. 

GR are usually associated with the first phonon of a 
small-amplitude harmonic motion. In the harmonic pic- 
ture, it corresponds to a coherent sum of one-particle one- 
hole (lplh) states [126]. GR are usually unbound and may 
decay by particle (p, n, d, a, 7...) emission, leading to an 
escape width in the GR spectra. The escape width is then 
due to a coupling of the correlated lplh states to the con- 
tinuum. Other contributions to the width of GR are the 
Landau damping and the spreading width. Landau damp- 
ing occurs due to a one-body coupling to non-coherent 
lplh states [134]. The spreading width is due to the resid- 
ual interaction coupling lplh states to 2p2h states. The 
2p2h states can also couple to 3p3h states and more com- 
plex npnh configurations until an equilibrated system is 
reached. These contributions to the GR width are illus- 
trated in Fig. 3.1. 

The proof of the vibrational nature of GR came with 
the observation of their two- and three-phonon states [135- 
137]. Multiphonon studies also provided a good test to 
the harmonic picture. In particular, anharmonicity was 
found in an abnormally large excitation probability of 
these states [135]. The coupling between different phonon 
states is predicted to be an important source of this an- 
harmonicity [138,139,60,140,61,64,141]. 

Coherent motion of fermions such as collective vibra- 
tions in nuclei can be modeled by time-dependent mean- 
field approaches like the TDHF theory. In fact, in its lin- 
earised version, TDHF is equivalent to the Random Phase 
Approximation (RPA) [142] which is the basic tool to 



Ar = Landau = Escape = Spreading 




Fig. 3.1. Illustration of the three contributions to the width 
of the GR (see text), (left) Landau damping due to coupling 
to incoherent lplh states, (middle) Escape width due to direct 
decay, (right) Spreading width due to coupling to 2p2h states. 



understand the collective vibrations in terms of indepen- 
dent phonons. In particular, time evolution of one-body 
(collective) observables, which can be estimated using a 
TDHF code, contains relevant informations on the vibra- 
tion properties, such as their energy spectra. 

Direct decay contributing to the escape width can be 
studied within the continuum- (Q) RPA model 3 [143-148]. 
GR direct decay can also be investigated within the TDHF 
framework [149-151,66]. Indeed, the TDHF theory is able 
to treat such process as it allows for the evaporation of un- 
bound components of the single particle wave-functions. 
Due to its one-body nature, Landau damping is also in- 
cluded in RPA and TDHF. However, they do not con- 
tain the residual interaction responsible for the spreading 
width. 

An interesting feature of TDHF applications is that 
they are not limited to small amplitude vibrations, unlike 
RPA, allowing for investigations of non linear effects in 
collective motions. In particular, couplings between collec- 
tive modes, which is a possible source of the anharmonicity 
discussed above, has been investigated [60,61,64]. 

In principle, TDHF and RPA codes can be used to 
study the vibrational spectra of any nucleus. However, 
most of the applications have focused on doubly magic 
nuclei. The main reason is that these nuclei are well de- 
scribed at the HF level, as there is usually no pairing cor- 
relation in their ground-state. Mid-shell nuclei, however, 
are better described with the HF+BCS or Hartree-Fock- 
Bogoliubov (HFB) approach. Recently, TDHF codes have 
been extended to study the role of pairing correlations on 
collective vibrations at the BCS level [114,66] and within 
the TDHFB theory. Indeed, TDHFB calculations with full 
Skyrmc EDF arc now possible in spherical symmetry [113], 
and more recently in three dimensions [65] . 

In this section, we first briefly introduce the linear re- 
sponse theory. For an illustrative purpose, applications of 
the latter to both low-lying vibrational states and GR are 
given within the TDHF framework. These standard ap- 
plications to the linear response theory are followed by 

3 The quasiparticle-RPA (QRPA) is an extension of the RPA 
including pairing correlations. 
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Name 


Approximation 


Variational 
space 


Associated 
observables 


TDHF 


mean-field (m.-f.) 


indep. part. 


one-body 


BV prescription 


m.-f. + small ampl. 
fluctuations 


indep. part 


one-body 
-(-fluctuations 


Path integrals 


m.-f. + fluctuations 


indep. part. 


one-body 
-^fluctuations 


(TD)GCM 


m.-f. + fluctuations 


correlated 
states 


one-body 
-(-fluctuations 


TDHF-Bogoliubov 


m.-f. + pairing 


indep. 
quasipart. 


generalised 
one-body 


Extended-TDHF 


m.-f. + collision 
(dissipation) 


correlated 
states 


one-body 


Stochastic- TDHF 


m.-f. + collision 
(dissipation-l- fluctuations) 


correlated 
states 


one-body 


Time-dependent 
density matrix 


m.-f. + two-body 
correlations 


correlated 
states 


one- and 
two-body 


Stochastic m.-f. 
(Functional integrals) 


Exact 

(within statistical errors) 


correlated 
states 


all 



Table 2.1. The TDHF approach and several possible extensions. 



an investigation of the direct decay by nucleon emission 
and its link to the GR microscopic structure. Then, non- 
linearities in collective vibrations arc studied to investigate 
a possible source of anharmonicity in GR multiphonon 
spectra. In the last application, we come back to the lin- 
ear response theory and study pairing vibrations with a 
TDHFB code. Finally we conclude this section and present 
some perspectives to the study of collective vibrations. 



3.2 Linear response theory 

The linear response theory has been widely used with 
TDHF to study collective vibrations in nuclei [104,152, 
153,149,150,154,60,57,58,56,155,63,156]. 

In this theory, one computes the time evolution of an 
observable Q(t) after an excitation induced by a small 
boost on the ground state 

|ff(0)) = e^^o) = |-?o)-|^^|^)+0(e 2 ), (3.1) 

where q v = (\P v \Q\)Po) is the transition amplitude between 
the ground state and the eigenstate \l? v ) of the Hamilto- 
nian with eigencnergy E v . The time-evolution of the state 
reads 

|#(t)) = e-^I^O)) 

= e- lE ^ h -jY^q^^W^ +0(e 2 ), 

(3.2) 

with hu) u = E u — Eq, 



The response Q{t) = (&(t)\Q\V(t)) - {& \Q\&o) to this 
excitation can be written 

Q{t) = -j^\q„\ 2 smut +()(<?). (3.3) 

The latter can be decomposed into various frequencies w, 
giving the strength function 

-H f°° 

Rn (w) = lim — / dtQit) sin(orf). (3.4) 
= J2\q v \ 2 6(u-u u ). (3.5) 

1/ 

We see in Eqs. (3.4) and (3.5) that the strength func- 
tion is obtained for small e. In practice, it is sufficient to 
check that the amplitude Q ma x of Q(t) evolves linearly 
with e [see Eq. (3.3)]. An example of evolution of Q m ax 
as a function of e is shown in Fig. 3.10-a in the case of a 
dipole response of the 132 Sn nucleus. 

It is interesting to note that the first phonon energy 
is obtained with an amplitude of the oscillation which is 
much smaller than the one associated to the first phonon. 
Indeed, in a coherent picture, and for a single mode with 
transition amplitude q, the number of excited phonons 
reads [60] 

(^Qmax \ 
— ) ' 

In particular, the amplitude associated to the first phonon, 
Qmax = 9/2 may be well beyond the linear regime. 

Fig. 3.2 gives an example for the octupole modes in 
208 Pb obtained with the tdhf3d code [55] (see also Ref. [57] 
for a study of this mode in 16 0). An octupole boost is ap- 
plied at initial time on the HF ground state of 208 Pb. The 
latter has spin-parity J v = + and the boost induces a 
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time (zs) 




20 

E (MeV) 

Fig. 3.2. (top) Time evolution of the octupole moment in 
208 Pb after an octupole boost obtained with the tdhf3d code 
for both TDHF (solid line) and unperturbed (dashed line) re- 
sponses (see text), (bottom) Associated strength function. 



transition with AL = 3, exciting vibrational states with 
J* = 3 _ . We can see on the upper panel of Fig. 3.2 the 
oscillation of the octupole moment induced by the boost 
(solid red line). The associated strength function, shown 
in the lower panel, exhibits a strong peak at an energy of 
~ 3.4 MeV. It corresponds to the main oscillation seen in 
the time evolution of Q30 (t) . Note that this state is clearly 
bound, as can be seen from the undamped nature of the 
oscillation. This peak is associated to the low-lying 3~ 
state in 208 Pb. The energy of this state is overestimated 
with the SLy4 parametrisation of the Skyrmc EDF, as the 
experimental value gives 2.6 MeV. However, its collective 
nature is unambiguous. This can be seen from a compar- 
ison with the unperturbed response of the same boost. 
The latter is obtained by freezing the mean-field in its ini- 
tial HF value, i.e., neglecting the self-consistency of the 
mean-field in the dynamics. This procedure removes the 
residual interaction which is responsible for the collectiv- 
ity of vibrations in TDHF (and RPA). We see that this 
peak disappears in the unperturbed spectrum, proving its 
collective nature. 

Another example is shown in Fig. 3.3. Here, a monopole 
boost is applied on the 208 Pb ground-state. The evolution 
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Fig. 3.3. (top) Time evolution of the mean square radius in 
208 Pb after a monopole boost obtained with the tdhfbrad 
code (without pairing), (middle) Time evolution of the oscilla- 
tion amplitude (log scale) . (bottom) Associated strength func- 
tion. From Ref. [157]. 



has been obtained with a TDHF code 4 in spherical sym- 
metry (allowing for large boxes) [157]. In this case, almost 
all the strength (if not all) goes into the GMR which is 
unbound. Such a resonance can then decay by particle 
emission. In particular, the direct decay induces an es- 
cape width. The latter is included in the TDHF framework 
thanks to evaporation of unbound components of the sin- 
gle particle wave functions [149-151]. Note that the emis- 
sion of nucleons in the continuum is quadratic with the 
boost intensity e. As a result, it does not affect the evo- 



4 This is the tdhfbrad code, developed by B. Avez, which 
solves the TDHFB equation. However, the pairing interaction 
is set to zero in the calculation shown in Fig. 3.3. 
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lution of Q(t) in the linear regime. However, special care 
should be taken with possible spurious numerical effects 
due to the reflection of the nucleons on the box edges [63] . 

The decay induces an exponential decrease of the oscil- 
lation amplitude (see middle panel of Fig. 3.3). The cor- 
responding escape width is = 0.8 MeV. This escape 
width is responsible for the main part of the width of the 
peak in the lower panel of Fig. 3.3. Indeed, the total width 
of the peak is r tot - ~ 1.1 MeV. The difference between the 
TDHF predictions of r tot - and T f could be attributed to 
Landau damping. 

Experimentally, the energy of the GMR in 208 Pb is [158] 

i^£*^ = 14.0 ±0.2 MeV, 
where the energy weighted moment mj is defined as [142] 




(3.6) 



10 15 
frequency co [MeV] 

Fig. 3.4. Dipole strength in the prolately deformed 238 U nu- 
cleus with the SkI3 parametrisation of the Skyrme functional. 
The strength is shown along the symmetry axis (z mode) 
and perpendicular to the symmetry axis (xy mode). The total 
strength is shown with the thin-solid line and compared with 
experimental data from Ref. [160] (thick-solid line). Adapted 
from Ref. [56]. 



This value is in excellent agreement with the TDHF result 
determined from the mi /mo ratio: E^f/^ F ~ 13.9 MeV [157] 
However, the total width obtained from TDHF, r^ HF ~ 
1.1 MeV, underestimates the experimental value -T*°*' = 
2.9±0.2 MeV [158]. This is due to the fact that TDHF does 
not take into account the spreading width (see Fig. 3.1 
and the discussion in the introduction of this section) . To 
go beyond and describe the spreading width, one would 
need to include the residual interaction. Calculations at 
the semi-classical level [159] and with the TDDM ap- 
proach [13, 12] have indeed shown an increase of the damp- 
ing thanks to the introduction of the collision term. 

TDHF calculations of giant-resonances are not limited 
to the study of monopole modes with spherical symmetry. 
Indeed, 3-dimcnsional TDHF codes have also been used 
to study other modes, such as the GDR and the GQR [60, 
57,56,58,64]. An example of such calculations is shown in 
Fig. 3.4, where the TDHF dipole strength function is plot- 
ted in 238 U for different excitation modes [56] and com- 
pared with experimental data from Ref. [160]. The 238 U 
being prolately deformed, the dipole response is investi- 
gated along and perpendicular to the deformation axis, 
exhibiting the well known splitting of the GDR peak, i.e., 
a lower (higher) energy along (perpendicular to) the de- 
formation axis [127]. The qualitative agreement between 
the theoretical prediction of the peak positions and exper- 
imental data is good. However, part of the TDHF width 
is due to the finite time window used in the Fourier trans- 
form. It is then difficult to draw any conclusion on the 
width of the GDR. 



«detector» 



ABC 



3.3 Direct decay of giant resonances 

Giant resonances usually lie above the proton and neu- 
tron emission thresholds. As mentioned above, their di- 
rect decay induces an escape width contributing to the 
total width of the GR. Such escape widths can be stud- 
ied with the TDHF formalism. Indeed, TDHF codes have 
been used in the past to investigate the direct decay of 



nucleus 

Fig. 3.5. Schematic description of the spatial repartition of nu- 
merical elements to compute spectra of emitted nucleons. The 
excited nucleus is in the center of a (spherical) box. Unbound 
parts of single-particle wave-functions are emitted in the con- 
tinuum. The "detector" shows the region of space where the 
energy of the emitted wave functions is computed. Absorbing 
boundary conditions (ABC) are used to absorb particles leav- 
ing the detector and to avoid spurious reflection on the box 
boundary. 



GMR [149,150]. In addition to its contribution to the es- 
cape width, the GR direct decay is particularly interesting 
as it brings informations on the microscopic structure of 
the GR [127]. This is the main purpose of the present 
section. 

TDHF calculations of GR direct decay are performed 
with large spatial grids to construct spectra of emitted nu- 
cleons with a good precision. It is then easier to use spher- 
ical TDHF codes, although the applications are limited to 
monopole vibrations only. Let us introduce a numerical 
" detector" corresponding to the region of space where the 
energy of the emitted nucleons is computed from Fourier 
transform of their spatial wave-functions. This detector 
should be away from the center of the box to avoid any 
nuclear interaction of the emitted nucleons with the nu- 
cleus. Fig. 3.5 shows a schematic representation of this 
numerical setup. Absorbing boundary conditions (ABC) 
with an imaginary potential may be used to avoid any 
spurious interaction with particles reflected on the box 
boundary [57,161]. 

Fig. 3.6 shows an example of a calculation of emit- 
ted nucleon spectra for the isovector GMR in 40 Ca. The 
calculations have been performed with the same spheri- 
cal code and simplified Skyrme EDF as in Ref. [162]. The 
upper panel of Fig. 3.6 shows the proton spectra at dif- 
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Fig. 3.6. Spectra of protons in the "detector" (see text and 
Fig. 3.5) at different times following an isovector monopole 
boost in 40 Ca. (top) Each spectrum is shifted vertically for 
clarity (time increases from bottom to top). The time delay 
between two consecutive spectra is AT = 5 zs. (bottom) The 
proton spectra at different times form an envelope (thick solid 
black line). The similar envelope obtained for neutron is also 
shown with a dashed line. 



ferent times. The first protons to reach the "detector" are 
obviously the fastest, i.e., with the highest kinetic energy. 
They also leave quickly the detector while slower protons 
reach it. Put together, these proton spectra form an en- 
velope 5 , shown in black thick solid line in the lower panel 
of Fig. 3.6. The neutron spectrum is also shown (dashed 
line). The fact that neutrons are more bound than protons 
in 40 Ca, together with the absence of Coulomb barrier for 
neutrons, explain that the neutron spectrum shows a more 
important contribution at low energy than the proton one. 

We can see in Fig. 3.6 that both proton and neutron 
spectra exhibit some structures which cannot be explained 
by a simple hydrodynamical model. Instead, one should 
seek for an explanation in terms of the microscopic struc- 
ture of the GR. This motivated a more detailed investi- 
gation with a realistic TDHF spherical code and a full 
Skyrme EDF [151]. A brief summary of the results for the 
GMR in 16 is presented here (see Ref. [151] for more 



5 This envelope is defined by the maxima of the spectra ob- 
tained at different times. 
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Fig. 3.7. (top) Time evolution of the monopole moment in 
16 O after a monopole boost obtained with the tdhfbrad 
code [113]. (bottom) Associated strength function and its de- 
composition onto single-particle quantum numbers I and j 
(spectroscopic notation). Solid (dashed) lines show neutron 
(proton) contributions. 



details and for more results on, e.g., tin isotopes). In this 
study, the tdhfbrad code [113] is used with the SLy4 
paramctrisation [49] of the Skyrme EDF without pairing. 
Fig. 3.7 shows the time evolution (top) of the monopole 
moment after an isoscalar monopole boost, and the as- 
sociated spectrum (bottom) obtained within the linear 
response theory. The GMR spectrum exhibits structures 
which are associated to different single-particle orbitals. 
For instance, the high energy shoulder around 31 MeV is 
due to Si/2 particle- hole excitations 6 . 

The spectra of emitted protons and neutrons arc shown 
in the upper panel of Fig. 3.8. The latter depend strongly 
on the associated single-particle quantum numbers. In par- 
ticular, no S1/2 nucleons are emitted. This is due to the 
fact that the lsi/2 hole state is deeply bound (-32.4 MeV 
for protons and -36.2 MeV for neutrons according to the 
HF initial configuration [151]). In fact, the high energy 
shoulder of the GMR spectrum (see Fig. 3.7-bottom) does 



6 The monopole excitations is associated to a AL — Oft an- 
gular momentum transfer so that particle and hole have the 
same quantum numbers at the time of the excitation. 
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Fig. 3.8. (top) Neutron (solid lines) and proton (dashed lines) 
direct-decay spectra, (bottom) Same spectra "shifted" by the 
energy of the initially occupied single-particle state. Their sum 
is shown in black solid line. 



not have enough energy to bring the initial lsi/2 particle 
into the continuum. 

The lower panel of Fig. 3.8 shows the same quantity 
as the upper panel, with a shift in energy (different for 
each lj contribution) corresponding to the binding energy 
of the hole state. The sum of each shifted lj contribution 
gives a spectrum which is very close to the GMR spectrum 
(compare with Fig. 3.7-bottom). The agreement is excel- 
lent for both the shape and the magnitude of the spectra. 
The origin of the structures in the direct emission spectra 
is entirely due to the shell structure of the nucleus. 

It is also interesting to note that, although the high 
energy shoulder in the 16 O GMR spectrum is due to the 
excitation of a bound lplh state, it appears in the "shifted" 
spectrum (Fig. 3.8-bottom) due to the emission of parti- 
cles inp-states. In fact, the TDHF (or RPA) residual inter- 
action is responsible for the coupling between the bound 
lplh si/2-state and unbound lpl/ip-states. Note that sim- 
ilar couplings have been obtained in tin isotopes [151]. For 
instance, the GMR in 100 Sn decays by protons only, while 
it is associated to a collective oscillation of both protons 
and neutrons. 

Coincidence experiments between particles emitted in 
the GR decay and the ejectile resulting from the GR exci- 
tation process have been performed in the past to investi- 
gate GR properties [127]. The present theoretical analysis 



of GR direct decay allows, in principle, a direct compar- 
ison between theoretical and experimental spectra. How- 
ever, for quantitative comparisons, one should use a more 
elaborated approach than the TDHF theory. Indeed, the 
fact that TDHF does not include lplh residual interac- 
tion is a strong limitation, as the latter has been shown 
to be crucial to reproduce the width and the fragmenta- 
tion of GR spectra [4] . The present analysis of GR decay 
should then be repeated with, e.g., the extended-TDHF 
(ETDHF) or the time-dependent density-matrix (TDDM) 
approaches (see table 2.1 and Ref. [3]). Note also that cal- 
culations should be performed with 2D or 3D codes in 
order to study the decay of GR with higher multipolarity 
than the GMR. 

3.4 Anharmonicity of giant resonances 

In the harmonic picture, a GR is the first phonon eigen- 
state of an harmonic oscillator describing the collective 
motion, and corresponds to a coherent sum of lplh states 
[126]. Experimental observations of 2 and 3-phonon states 
proves the vibrational nature of GR. However, they also 
show limitations of the harmonic picture [135-137]. In par- 
ticular, the excitation probability of multi-phonon states 
is larger than predicted by the harmonic picture. This in- 
dicates that different phonon states may be coupled by 
the residual interaction [138,139,60,140,61,141,64]. 

The TDHF approach takes into account some effects of 
the residual interaction if the considered phenomenon can 
be observed in the time evolution of a one-body observ- 
able. In particular, the nonlinear response in TDHF con- 
tains the couplings between one- and two-phonon states 
coming from the iplh and Ipih residual interaction [60]. 
In that sense, it goes beyond the RPA, which is a harmonic 
picture and contains only lplh residual interaction. 

The couplings leading to the excitation of a GMR or 
a GQR (resp. a GMR) on top of a GDR (resp. a GQR) 
were investigated in Ref. [60] using the nonlinear response 
to an external field in the TDHF theory. As a continuation 
to this work, different techniques to compute the matrix 
elements of the residual interaction responsible for these 
couplings were introduced in Ref. [64]. 

These matrix elements can be written = (z^l/lz^) 
where the residual interaction V = H — Hq is the differ- 
ence between the full Hamiltonian H and the HF+RPA 
Hamiltonian Hq. \v) and are 1 and 2-phonon eigen- 
states of Hq with eigenenergies E v = Eq + hjj v and E v ^ t = 
Eq + hujjy + huj^ 1 where oj„ iM denote the collective frequen- 
cies and Eq is the ground state energy. The state can 
be seen as one phonon of the GR /1 (e.g., a GQR) excited 
on top of one phonon of the GR v (e.g., a GDR). 

In addition to the original technique based on the non- 
linearities of the time-dependent response [60], two other 
methods were introduced in Ref. [64] . A brief summary of 
these three methods is given below (see Ref. [64] for more 
details): 

— method 1: A boost eT applied on the ground state 
induces, at lowest order in k v , an oscillation of Q v {t) = 
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Fig. 3.9. Time evolution of the dipole (a) and quadrupole (b) 
moments in 132 Sn after a dipole boost. 



{Qv){i) linear in k v , and an oscillation of quadratic 
in k v and proportional to v^. Computing the response 
Qn(t) to such a boost with TDHF gives then access 

to Vft. 

method 2: The same boost can be applied on a HF 
state obtained with a small constraint XQ^. The linear 
response Q v {t) oscillates then with a frequency w„(A). 
The variation is proportional to v^. Computing 
^(A = 0) with a TDHF or a deformed RPA code 7 
gives also access to v^. 

method 3: A "double" boost e~ lk "®" e~ lk vQv \ s applied 
on the ground state |0). Define the "coupling response 
function" as 8 



0.1 



0.2 

k D (fm" 1 ) 



Fig. 3.10. Circles indicate the first minimum and maximum of 
the dipole (a) and quadrupole (b) moment evolution, respec- 
tively, following a dipole boost in 132 Sn with a boost veloc- 
ity ko- Dashed lines show linear and quadratic extrapolations 
at ko — > of the dipole (a) and quadrupole (b) amplitudes, 
respectively. 



amplitude of the dipole (resp. quadrupole) oscillation is 
linear (quadratic) in the boost velocity kjj. Numerical 
application gives a matrix element of the residual inter- 
action t)!!' 



-0.61 MeV. The two other methods give 



,(2) 



-0.56 and v, 



(3) 

Q 



-0.68 MeV, respectively, show- 



-1 



dt cos(u)t)Q v (t). 



(3.7) 



At lowest order in k^ v , we can show that is 
proportional to and exhibits peaks at u> u and \w v ± 
Wu\ with opposite signs. The amplitude of these peaks 
provides a third way to extract v^. 

An illustration of the first method applied to the 132 Sn 
nucleus with v =GDR and \i =GQR is given in Figs. 3.9 
and 3.10. On Fig. 3.9, we observe an oscillation of both 
the dipole and quadrupole moments, although the boost 
contains only the dipole moment. The oscillation of the 
quadrupole moment is, in fact, induced by the residual 
interaction. The right panel shows that, as expected, the 



7 This technique involves linear response only. The matrix 
element can then be computed with a RPA code allowing 
initial deformations (generated by the constraint XQ) of the 
vacuum. 

8 This function differs from the standard response function 
[see Eq. (3.4)] essentially by the cosine function instead of a 
sine function. 



ing a relatively good agreement between the three meth- 
ods [64]. 

Couplings have been computed in other tin isotopes [64] 
and in other nuclei ( 40 Ca, 90 Zr, and 208 Pb) with TDHF [60] 
and with a boson mapping method [140]. Refs. [60,140] 
also discuss couplings involving the GMR built on top of 
the GQR or the GDR. The TDHF results provide a confir- 
mation to the amplitude of the couplings computed with 
the boson mapping method (see discussion in Ref. [60]). 

Another conclusion of Ref. [64] is that there is no 
(or little) dependence of the coupling between dipole and 
quadrupole motion with isospin. However, an overall de- 
crease of the coupling is obtained with increasing mass, 
indicating that the couplings are mediated by the sur- 
face [60,140,64]. 

Overall, the couplings are small but significant com- 
pared to the GR energies (e.g., v/u ~ 5% for 132 Sn). 
Their effect on the first phonon is negligible, but becomes 
sensible on the second and third phonon, with a typical 
shift in hid of the order of ~ 0.5 MeV as compared to 
the harmonic picture [140]. How the anharmonicitics in- 
duced by these couplings affect the excitation probabil- 
ity to the multiphonon states have been investigated by 
Lanza et al. within a semiclassical coupled-channels for- 
malism [141]. This model, based on the boson mapping 
method for the multi-phonon properties, allows for cal- 
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culations of inelastic cross sections for the multiple ex- 
citation of giant resonances induced by heavy-ion probes. 
Their calculations show that these anharmonicities induce 
an increase of the inelastic cross-section (as compared to 
the harmonic model) in the multi-phonon region, in good 
agreement with experimental data. 

The role of pairing correlations, neglected so far, should 
also be considered. For instance, fully self-consistent quasi- 
particle-RPA (QRPA) codes allowing for static deforma- 
tions (see, e.g., the code developped by S. Peru [163]) could 
be used to obtain the couplings between the GQR and the 
GMR built on top of it. Couplings with exotic modes such 
as the pygmy dipolc resonance [164] should also be inves- 
tigated with the present methods. 



3.5 Pairing dynamics 

The TDHF calculations presented in the previous sections 
were applied to "normal" vibrations, i.e., vibrations of the 
one-body density (also called normal density) p(t) with 
matrix elements 




p a p{t) = (#(t)|aja„|#(t)>. 



(3.8) 



These vibrations do not probe directly the pairing correla- 
tions between nucleons. Inclusion of pairing is possible in 
the small amplitude limit with the QRPA based on HFB 
vacua. The HFB +QRPA has been widely used in nuclear 
structure studies [165,148,166,163]. 

Similarly to the fact that the TDHF approach is an 
extension to the HF+RPA, a natural extension to the 

HFB+QRPA is the time-dependent Hartree-Fock-Bogoliubov^hanks to a gauge symmetry breaking. The Hartree-Fock- 



2p 



Fig. 3.11. Illustration of the configurations used to generate a 
ground state with pairing correlations, (top) Exact case, (bot- 
tom) BCS approximation. 

mostly (but not only) time-reversed states. The pairing 
residual interaction induces a scattering of a pair of nucle- 
ons across the Fermi surface. As a result, the ground-state 
with pairing correlations is a sum of 2p2h states where the 
2p (2h) are essentially time-reversed states. Such a state is 
represented schematically in the upper part of Fig. 3.11. 

The treatment of pairing correlations in finite nuclei 
is simplified with the Bardeen-Cooper-Schrieffer (BCS) 
approximation initially developed to interpret supracon- 
ductivity in metals [169]. In this approximation, the 2p2h 
states are replaced by a sum of 2p and 2h states (see bot- 
tom part of Fig. 3.11). The resulting approximation of the 
ground-state can then be written as a vacuum of quasi- 
particles, allowing for the application of the Wick theo- 
rem [142]. As a result, the BCS approximation leads to a 
generalised mean-field theory. 

As we can see in Fig. 3.11, the "price to pay" is that the 
BCS ground-state is not an eigenstate of the particle num- 
ber operator anymore. Pairing correlations are included 



(TDHFB) theory [2]. In particular, the TDHFB theory 
provides a fully self-consistent 9 response to an external 
excitation including pairing dynamics and non-linearities. 

The development of a realistic TDHFB spherical code 
with a full Skyrme EDF and a density-dependent pair- 
ing effective interaction was initiated in Rcf. [113] and 
applied to the study of pairing vibrations. Recently, 3- 
dimcnsional codes have been developed to study the effect 
of the pairing interaction on "normal" vibrations at the 
BCS level [114,66] and solving the TDHFB equation with 
a Skyrme functional [65] and with the Gogny effective in- 
teraction [48]. The present section introduces briefly the 
TDHFB formalism and discusses the application to pair- 
ing vibrations presented in Ref. [113]. 



3.5.1 The TDHFB theory 

Pairing correlations are essentially due to an attractive 
short range contribution of the residual interaction in the 
1 Sq channel 10 [167,168]. Pairing correlations affect then 



Bogoliubov (HFB) theory is more general but shares the 
same features with the BCS approach. In fact, it is an ex- 
tension to the BCS theory where pairs are not limited to 
time-reversed states. 

In the (TD)HFB theory, all the information on the 
state of the system is contained in the generalised density 
matrix 1Z defined as 



K 



((S 3 -Oi)) ' 



P 

— K* 



1 



(3.9) 



where k is the so-called pairing tensor, k and n* contain 
the pairing correlations (at the HFB level). 

The time evolution of the generalised density matrix 
is given by the TDHFB equation 

ih^ = [U,K] (3.10) 

which has the same form than the TDHF equation (2.5). 
The generalised HFB Hamiltonian reads 



9 Here, the self-consistency refers to the fact that the HFB 
vacuum and the residual interaction inducing the collective 
dynamics are derived from the same EDF. 
10 This notation means that the two nucleons are coupled to 
produce a total isospin 1, a total orbital angular momentum 
L = (S— wave), and a total spin 0. 



where 



H = 



S£[p, K, K* 
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-A* 
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-h* 
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SS[p, 



(3.11) 
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are the HF Hamiltonian and the pairing field, respectively, 
and £ [p, «, k*] is the EDF including pairing. 



3.5.2 Application to pairing vibrations 

Pairing vibrations are a particular manifestation of the 
dynamics of pairing correlations [126,142,170]. They are 
probed in two-nucleon transfer reactions [171-175]. Pair- 
ing correlations are then expected to induce a collectivity 
which manifests itself as an increase of transition ampli- 
tude toward states associated to pairing vibrations. Start- 
ing with an even-even nucleus ground-state with A nucle- 
ons and spin-parity + , and assuming a AL = direct 
pair transfer reaction, pair vibration states with J 71 = + 
are populated in the A + 2 (pair addition) and/or A — 2 
(pair removal) nuclei. 

Such a process can be simulated within the TDHFB 
formalism using an initial boost with a Hcrmitean pair- 
transfer operator [170] 



F= / dr/(r) 



'r.t 



+ a rj -|-a r ,^ 



(3.13) 




600 
time (fm/c) 
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where the arrows label the spin of the single-particles (we 
omit the isospin to simplify the notation). In the present 
application, f(r) is a Fermi-Dirac spatial distribution con- 
taining the nucleus and cutting at 4 fm outside the nu- 
cleus. Its role is to remove unphysical high energy modes 
associated to pair creation outside of the nucleus. 

In this approach, it is assumed that the spectroscopy of 
the A— 2, A and A+2 nuclei can be obtained from the same 
quasiparticle vacuum (the A ground-state). Note that re- 
cent improvements have been proposed by Grasso et al. 
where this limitation is overcome for ground-state to ground- 
state transitions by using different vacua for the parent 
and daughter nuclei [176]. 

The tdhfbrad code has been developed to solve the 
TDHFB equation in spherical symmetry with a full Skyrme 
EDF and density-dependent pairing effective interaction. 
As a first application, the linear response of (F)(t) has 
been computed in several oxygen and calcium isotopes [113] . 
The timc-cvolution of (F) (t) is shown in the upper panel 
of Fig. 3.12 for a 18 vacuum. The apparent chaotic be- 
haviour of (F) (t) is due to the simultaneous excitation of 
several pair vibrations, as we can see from the strength 
function (solid line) in the middle panel. Both pair addi- 
tional and pair removal (indicated by the arrows) modes 
are present. A comparison with the unperturbed strength 
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Fig. 3.12. (top) Evolution of (F)(t) after a pair transfer type 
excitation on ls O. The inset shows the same quantity at early 
times, (middle) Associated TDHFB strength function (solid) 
compared with the unperturbed spectrum (dashed). The ar- 
rows indicate pair removal transitions from the lps/2 (solid) 
and lpi/2 (dotted) deep hole states, (bottom) TDHFB strength 
function decomposed into single-particle orbital angular mo- 
mentum i-components. 



The bottom panel of Fig. 3.12 shows a decomposition 
of the response in terms of single-particle orbital angu- 
lar momentum I. Together with the structure of the ini- 
tial HFB vacuum, this decomposition allows for an un- 
derstanding of each peak in terms of their main particle 
and hole contributions. See Rcf. [113] for a detailed mi- 



function (dashed line) obtained by removing the self-consistency ic ^ ^ gee algQ Refg " [n3 m] for an ^ ig 



of the generalised mean-field shows two features in the 
TDHFB spectrum: 

— an increase of the strength, 

— and a lowering of the transition energies. 

Both are compatible with the attractive nature of the dy- 
namical pairing residual interaction. In particular, the in- 
crease of the strength is a clear signature for collective ef- 
fects. Note that similar conclusions were drawn from con- 
tinuum QRPA calculations by Khan et al. [173]. 



of other nuclei (oxygen, calcium and tin isotopes). 

This first realistic application of the TDHFB theory 
to nuclear systems has confirmed previous QRPA calcula- 
tions of pairing vibrations [173]. Future applications of the 
tdhfbrad code to study the decay of GMR, and to in- 
vestigate non-linear effects in pairing dynamics are envis- 
aged. The recent development of 3-dimcnsional codes [65, 
114,48] opens also new perspectives for the study of L ^ 
vibrations. 
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3.6 Conclusions and perspectives 

Real time mean-field calculations have been performed to 
investigate collective vibration properties. The response 
to an external excitation has been obtained with 3D and 
spherical TDHF codes and associated strength functions 
have been computed within the linear response theory. 

Direct decay of GR have been analysed from energy 
spectra of emitted nucleons. Within the TDHF approach, 
the latter contains enough information to reconstruct the 
strength function if the hole structure of the nucleus is 
known. A comparison between the microscopic decompo- 
sitions of the strengths obtained from the time response 
of the excitation operator and from the emitted nucleon 
spectra shows that the residual interaction couples bound 
lplh states with unbound ones, allowing for particle emis- 
sion in the continuum. 

Non linear vibrations were also studied within the TDHF 
framework. They are used to quantify the coupling be- 
tween one GR phonon and two (different) GR phonon 
states. The large values of the couplings which have been 
obtained in different nuclei confirm that these kind of cou- 
plings is a probable source of the anharmonicitics observed 
in GR multiphonon experiments. 

Recents works were devoted to the inclusion of pairing 
correlations in the mean-field dynamics. To that purpose, 
a realistic spherical TDHFB code has been built and ap- 
plied to the study of pairing vibrations excited in pair- 
transfer reactions. Comparisons with unperturbed calcu- 
lations show that the dynamical pairing residual interac- 
tion included in TDHFB is attractive and induces some 
collectivity to the pairing vibrations. 

Several possible extensions to these works have already 
been discussed, such as the study of non-linear vibrations 
and particle decay with pairing, and the study of L ^ 
vibrations with 3-dimensional TDHFB codes. 

Other perspectives could be considered. 

In the present studies, we focused on collective mo- 
tion at zero temperature, in particular vibrations built on 
top of the ground state. The role of finite temperature on 
collective motion (e.g., the so-called "hot GR") has been 
widely discussed in the nuclear physics community [178]. 
Questions such as the effect of temperature on pairing 
dynamics and on the couplings between GR multiphonon 
states could be addressed with an extension of the present 
calculations to finite temperature systems. TDHF studies 
of giant resonances at finite temperature are indeed possi- 
ble [179] starting from an initial hot HF solution [180]. In 
section 4.2.3, we also discuss the particular case of GDR 
excited in the fusion of two nuclei. 

Recently, it has been shown with a molecular dynamics 
approach that GR in light nuclei ( 12 C, 16 O, 24 Mg) were 
also affected by a-clustering [181]. In particular, these a- 
clustcrs induce new vibrational modes which could cou- 
ple to "standard" GR and produce a new source of an- 
harmonicity. The coupling between these new modes and 
standard vibrations could be investigated in the future 
with the present techniques. 

The calculations presented in this section and the pos- 
sible perspectives discussed above are based on a mean- 



field approach. For a more realistic comparison to experi- 
mental data, extensions to theories going beyond the one- 
body limit are mandatory, in particular to reproduce the 
GR fragmentations and widths [4] . Possible approaches in- 
clude extended-TDHF [182], second RPA [183,184], time- 
dependent density matrix [14-16,185], or stochastic one- 
body transport [121,123,125] theories. Although they all 
face their own technical difficulties, these approaches could 
benefit from the recent increase of computational power. 



4 Heavy-ion collisions 

The previous section dealt with small amplitude vibra- 
tions. We now investigate the collision of two atomic nuclei 
which is by nature a large amplitude collective motion. We 
first start in section 4.1 by describing some motivations to 
study heavy-ion collisions. 

The outcome of a heavy-ion collision depends essen- 
tially on few properties of the entrance channel: energy, 
masses, angular momentum, deformation and orientation, 
N/Z asymmetry, and internal structure (e.g., magicity). 
Thus, we discuss different reaction mechanisms and their 
dependence on entrance channel properties. In sections 4.2, 
4.3, 4.4, and 4.5, we present studies on fusion, transfer, 
deep-inelastic, and quasi-fission reactions, respectively. In 
section 4.6 we present a theoretical study of actinidc col- 
lisions. Note that only a brief summary of each study is 
given, in particular when more details can be found in 
the associated publications. Finally, we conclude and give 
some perspectives. 



4.1 Motivations 

There are several motivations to achieve a good under- 
standing of heavy-ion reaction mechanisms. Two of them 
are particularly relevant to the work presented in this sec- 
tion: 

— To investigate quantum phenomena with complex sys- 
tems, 

— To optimize the production of specific nuclei. 
Some examples are also discussed below. 



4.1.1 Unique perspectives brought by the nuclear case to 
study quantum phenomena 

A huge difference between the quantum treatment of a 
macromolecule, such as a fullcrene, and an atomic nucleus 
is their interaction with their environments. Indeed, large 
molecules interact strongly with the surrounding gas and 
photon bath, while a system of colliding nuclei is free of 
such interactions. This is because (i) nuclei do not inter- 
act with most photons in nature due to their small size 
(few fm), and (ii) nuclear excited states have much longer 
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lifetimes (~ 10 -18 s) than typical reaction times in heavy- 
ion collisions 11 (< 10~ 20 s). Being entirely isolated during 
the reaction, heavy-ion collisions could then be considered 
as ideal to investigate the role of the composite nature of 
nuclei on quantum processes such as tunnelling 12 . 

As an example, one of the most striking signatures 
of quantum phenomena in nuclear collisions is the coher- 
ent quantum superposition of reaction channels leading 
to fusion barrier distributions [187] which are interpreted 
within the coupled-channel framework [188]. Fusion via 
tunnelling is effectively a powerful microscope, magnifying 
quantum effects in nuclear collisions through the exponen- 
tial dependence of tunnelling probability on the potential, 
which is modified by quantum coherent effects. 

4.1.2 Production of specific nuclear systems: case of the 
heaviest nuclei 

Stellar nucleosynthesis has produced a great variety of nu- 
clei thanks to various reaction mechanisms. Of course all 
that remains on Earth are the stable and very long-lived 
by-products of these reactions, and in order to perform ex- 
perimental nuclear studies, these must serve as our start- 
ing points. The internal structure of nuclei is of consider- 
able interest in their own right and we can use the nuclei 
left by nature to create others, like new super-heavy ele- 
ments (SHE) or isotopes which may or not be produced 
at the far end of the astrophysical r-process [189]. 

Transfermium nuclei are of great interests for our un- 
derstanding of the quantum nature of the nuclear many- 
body problem as most of them are stable against fission 
thanks to quantum shell effects only. In particular, an is- 
land of stability is predicted in the super-heavy region 
with Z = 114-126 and N = 172-184 (see, e.g., Ref. [190] 
for mean-field calculations). 

The discovery of SHE is important for various funda- 
mental research fields: 

— Astrophysics: Does the r-process reach the super-heavy 
island of stability, or does it stop with spontaneous 
fissile lighter nuclei (see, e.g., [189])? 

— Atomic physics and chemistry: Relativistic effects on 
the electrons play a role on the atomic structure of 
SHE. As a result, SHE chemical properties may devi- 
ate from the standard classification suggested by the 
Mendeleiev table (see, e.g., [191,192]). 

— Nuclear physics: The precise location of the island of 
stability at the top of the nuclear chart would constrain 
the theoretical models (see, e.g., [190]). 

The most important breakthrough in the quest for 
SHE came from heavy-ion fusion-evaporation reactions [193 

11 This point is crucial. Indeed, if an excited state decays by 
7-emission during the collision, then the state of the environ- 
ment changes from the vacuum 07-state to 17-state. If such an 
emission occur, the system of nucleons would then be strongly 
coupled to its environment. 

12 If the system was not isolated from its environment, then 
decoherence would occur and the system would have a classical 
behaviour [186]. 



Elements up to Z=118 have been produced with this tech- 
nique [194]. However, the cross-sections for the production 
of SHE in their ground state are extremely small. The 
smallest cross-section measured is 30 femtobarns (fb) with 
the production of two Z = 113 nuclei, after more than 
half-year irradiation of 209 Bi targets with 70 Zn beam [195]. 

Two main factors are hindering SHE formation by fu- 
sion: 

— The quasi-fission mechanism strongly dominates [196- 
198]. It leads to a fast re-separation of the projectile 
and target-like fragments after an exchange of nucleons 
from the heavy fragment to the light one. 

— In case of fusion and formation of a compound nucleus 
at finite energy, the survival probability against fission 
is very small due to the small fission-barriers of SHE. 

A summary of the heaviest nuclei, most of them pro- 
duced artificially on Earth, is given in Fig. 4.1. All trans- 
fcrmiums formed by fusion of stable nuclei are neutron- 
deficient. As a consequence, nothing is known about j3- 
stable or neutron-rich transfermiums. 

There are three possible mechanisms which could pro- 
duce neutron-rich transfermium nuclei and, then, get closer 
to the island of stability: (i) a rapid neutron capture pro- 
cess, (ii) fusion with neutron-rich nuclei, and (Hi) multi- 
nucleon transfer in actinide collisions. Neutron fluxes in 
nuclear reactors are not sufficient to use (i). A critical lim- 
itation of the upcoming exotic beam facilities is the beam 
intensities, without which the feasibility of (ii) is ques- 
tionable. However, the process (Hi), which is discussed in 
Sec. 4.6, seems promising. The main advantages of using 
actinides are: 

— They exist in nature (e.g., 238 U) or can be produced 
as radioactive targets (e.g., 248 Cm). 

— They have ~ 50% more neutrons than protons, which 
is of great help for the production of new isotopes along 
the /3-stability line. 

4.2 Fusion with medium mass systems 

Dynamical mean-field theories like TDHF are well suited 
to the study of low-energy reaction mechanisms, such as 
fusion, at energies around the barrier. Indeed, at low ener- 
gies, the Pauli principle blocks collisions between nucleons, 
increasing their mean-free path to the order of the size of 
the nuclear system. In addition, fusion occurs by transfer- 
ring relative motion into internal excitation via one-body 
mechanisms well treated by the TDHF approach. 

Early TDHF codes have been successfully applied to 
describe above-barrier fusion reactions in light systems [200] 
' However, these calculations also predicted a lower limit 
to the angular momentum for fusion. For smaller angular 
momenta, a "transparency" was observed in the calcula- 
tions (two fragments are emitted along the collision axis) . 
This prediction was never confirmed experimentally. In 
fact, it was shown by Umar and collaborators that this 
so-called " fusion-window" problem was solved with the in- 
clusion of the spin-orbit interaction [153]. Indeed, the lat- 
ter was shown to be an important source of dissipation in 
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Number of neutrons N 

Fig. 4.1. Upper part of the nuclear chart. The dark blue area indicates quantum shell correction energy leading to more 
bound, and then more stable nuclei, according to a microscopic-macroscopic approach [199]. The purple line shows the expected 
/^-stability line. The expected path followed by the r-process is shown with a red line. 



heavy-ion collisions. Modern TDHF calculations are now 
performed with a full Skyrme EDF including spin-orbit 
terms [55-58]. 

The first observable we propose to study is the cap- 
ture threshold between two nuclei. We consider here sys- 
tems with Z1Z2 < 1600 to avoid the well known fusion 
hindrance observed in very heavy systems 13 . We also dis- 
cuss above barrier fusion cross-sections, and investigate 
the path to fusion from the charge equilibration process. 



4.2.1 Fusion of spherical nuclei: the example of the 
i6 + 208p b system 

A reference nucleus- nucleus potential could be obtained 
in the frozen approximation with HF (or HFB) densi- 
ties [201], where the collision partners arc assumed to 
keep their ground-state density during the approach. The 
frozen potential can be computed with the same Skyrme 
EDF as in the TDHF calculations by translating the nu- 
clei in their HF state [202]. Comparisons between TDHF 
and frozen fusion barriers allow to identify the role of dy- 



1 The case of heavier systems exhibiting fusion hindrance is 
discussed in section 4.5. 



namical effects, which are included in TDHF but absent 
from the frozen approach. 

Writing the HF energy E[p] as an integral of an energy 
density 'H[/ci(r)], i.e., 

E[p] = J dr H[p(r)}, (4.1) 

we get the expression for the frozen potential 

V(R)= fdr H[pi(T)+p 2 (r-Ii)}-E[p 1 \-E[p 2 }, (4.2) 

where R is the distance between the centers of mass of the 
nuclei, and pi.2 are the densities of their HF ground-state. 
In order to take into account the Pauli principle between 
the nucleons of one nucleus and the nucleons of the other 
one in Eq. (4.2), a proper treatment of the kinetic energy 
has to be considered, for instance using the Thomas- Fermi 
approach [203,204]. However, for light and medium heavy 
systems, the barrier radius is large enough to neglect the 
Pauli principle between the two reactants at the barrier. 
We then neglect the latter in the determination of the 
HF- frozen barrier height. 

An example is shown in Fig. 4.2 for the 16 O+ 208 Pb 
system which could be considered as a benchmark in low 
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Fig. 4.2. Nucleus-nucleus potentials obtained with the frozen 
approximation in the 16 O+ 208 Pb system. Nuclear (dashed blue 
line) and Coulomb (dotted red line) contributions, and their 
sum (solid line) have been obtained with the frozen approxi- 
mation. For comparison, the Wong potential [208] is shown in 
green dot-dashed line. The latter is obtained with a potential 
depth Vq = 70 MeV, a potential diffuseness a = 0.48 fm, and 
nuclear radii Ri — l.25A]^ 3 fm. 



energy reaction studies (see, e.g., Refs. [205-207,105,100, 
99]). A comparison with the Wong formula [208] is given. 
We see that, although the potential heights agree, differ- 
ences appear at short distances. In fact, inside the barrier, 
the parametrisations of nucleus-nucleus potentials are less 
constrained by experimental data. In addition, as men- 
tioned before, the frozen approximation neglects the Pauli 
principle for the nucleons of different collision partners, 
which may affect the inner barrier region. 

The height of the barrier obtained with the frozen 



approximation is V, 



r frozen 
B 



76.0 MeV. This value is 
close to the barrier obtained with the Wong formula [208] , 
yWong _ ?5 g MeV ^ whUe it ig 1 MeV smaller than the 

Bass barrier [210], V^ ass ~ 77.0 MeV. All these barri- 
ers overestimate the experimental value obtained from the 
centroid of the barrier distribution (see Fig. 4.3) Vg Xp ' ~ 
74.5 MeV. 

To investigate the possible role of dynamical effects 
on the fusion barrier, the latter has been computed with 
the tdhf3d code [55]. Because the TDHF theory does 
not allow for quantum tunnelling of the many-body wave 
function, the TDHF fusion barrier can be identified as 
the capture threshold for central collisions, above which 
a compound system is formed and below which the two 
fragments re-separate. Due to the finite time of the TDHF 
evolutions, one has to define a maximum computational 
time 15 above which the final configuration (i.e., one com- 
pound system or two fragments) is assumed to be reached. 



14 The same value has been obtained, independently, by 
Washiyama et al. [78] and Guo et al. [209]. 

15 This time may depend on the system. For medium mass 
systems such as 16 O+ 208 Pb, a typical time of 10 3 fm/c 
(1 zs=300 fm/c) is used. 
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Fig. 4.3. Experimental fusion barrier distribution of the 
i6 O+ 208p b system from Ref [207]. The frozen, Bass and 
TDHF barriers (see text) are shown with arrows. 
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Fig. 4.4. Relative distance between the fragments as a func- 
tion of time for head-on 16 O+ 20S Pb reactions computed with 
the tdhf3d code [55]. 



Figure 4.4 shows the evolution of the relative distance 
between fragment centers of mass in central 16 O+ 208 Pb 
collisions at different energies around the capture thresh- 
old. We clearly identify two sets of trajectories associated 
to capture (fusion) and to re-separation of the fragments. 
These calculations predict a fusion threshold of VfJ DHF = 
74.445 ± 0.005 MeV. As a result, the dynamical effects 
included in TDHF calculations lower the barrier by ~ 
1.5 MeV for this system as compared to the HF- frozen cal- 
culation. We observe in Fig. 4.3 a good agreement between 
the TDHF prediction and the centroid of the experimen- 
tal barrier distribution. Other methods based on a macro- 
scopic reduction of the mean-field dynamics, namely the 
dissipative-dynamics TDHF [78] and the density-constrained 
TDHF [80], also find similar results with an energy de- 
pendence to the barrier heights ranging from 74.5 at low 
energies to 76 MeV at higher energies where the frozen 
approach is expected to be more reliable. 
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Fig. 4.5. (top) Density evolution for the reaction 16 O+ 208 Pb 
corresponding to a head-on collision at a center of mass energy 
Ec.m. = 74.44 MeV (just below the fusion barrier). The red 
surfaces correspond to an iso-density at half the saturation 
density (p /2 = 0.08 fm~ 3 ). Each fi gure is separated by a time 
step of 135 fm/c. Time runs from left to right, (bottom) Same 
at E c , m , = 74.45 MeV, i.e., just above the fusion threshold. 
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Fig. 4.6. Bass barriers [210] (stars) and barriers extracted 
from TDHF calculations (circles) as a function of experimental 
barriers (centroids of fusion barrier distributions). 
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Fig. 4.7. Experimental fusion cross sections from Ref. [207] 
(stars) compared to cross sections deduced from TDHF cal- 
culations (lines) of 16 O+ 208 Pb collisions. The two lines cor- 
respond, respectively, to lower and upper limits of theoretical 
cross sections. 



has been obtained (better than the Bass parametrisation) 
for all the studied systems. Other calculations with 3- 
dimensional TDHF codes confirmed the predictive power 
of the TDHF approach for the determination of fusion 
barriers [78,209]. 

Above barrier fusion cross-sections have been computed 
for the 16 O+ 208 Pb system in Ref. [76]. The fact that fu- 
sion probabilities are either or 1 implies that cross sec- 
tions are obtained using the "quantum sharp cutoff for- 
mula" [211] 



To get a deeper insight into these dynamical effects, 
the density evolutions at E c , m , = 74.44 and 74.45 MeV 
are plotted in Fig. 4.5. In the first case, a "di-nuclear" sys- 
tem is formed during a relatively long time (~ 500 fm/c) 
before re-separation. In the second case, the system over- 
comes the fusion barrier. More generally, the two figures 
illustrate the richness of physical phenomena contained in 
TDHF calculations: surface diffuseness, neck formation, 
quadrupolc/octupole shapes of the compound system... 

The observed lowering of the fusion barrier due to dy- 
namical effects could be partly explained by a coupling of 
the relative motion to a transfer mechanism [76]. In fact, 
the outgoing channel of 16 O+ 208 Pb at E c . m . = 74.44 MeV 
(see top of Fig. 4.5) is, in average, 14 C+ 210 Po. This two- 
proton transfer channel effectively lowers the barrier by 
decreasing Z1Z2 and, then, the Coulomb repulsion. Trans- 
fer reactions in the 16 O+ 208 Pb system are discussed in 
more details in section 4.3. Note that low-lying collective 
vibrations, such as the first 3~ state in 208 Pb (see Fig. 3.2) 
also affect the fusion barrier distribution [207]. 

Systematic calculations of fusion barriers have been 
performed for medium mass systems involving spherical 
nuclei [76]. A summary of the results is shown in Fig. 4.6. 
A good reproduction of the barrier distribution centroids 



cT fus (E) = ^[l max (E) + l} 2 , (4.3) 

where the fusion probability is for I > l ma ,x(E) and 1 
for I < l m ax{E). To avoid discontinuities due to the in- 
teger values of l max {E), [l max (E) + l]h is gener ally a p- 
proximated by its semi-classical equivalent C c = \J2jiE b c . 
The latter corresponds to the classical angular momentum 
threshold for fusion and b c denotes the maximum impact 
parameter below which fusion takes place [212]. We fi- 
nally obtain the standard classical expression for fusion 
cross sections <jf us (E) ~ TrC^/2fiE = 7t6 2 . 

The results are shown in Fig. 4.7 for the 16 O+ 208 Pb 
system. Fusion cross-sections are overestimated by about 
16% above the barrier. Although this discrepancy is small 
for a theory which has no parameter adjusted on reaction 
mechanisms, its origin is unclear. 

Finally, the calculations are not able to reproduce the 
sub-barrier energies. This is of course one of the main 
drawbacks of the TDHF approach. The inclusion of quan- 
tum tunnelling of the many-body wave function is clearly 
one of the biggest challenges in the microscopic treatment 
of low-energy nuclear reactions. 
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Fig. 4.8. (left) Experimental barrier distributions for 

16 + 154g m j 213 ] ( rig]lt ) game for the 16 + 238 U system [ 2 14]. 

The arrows indicate the barriers obtained from TDHF calcu- 
lations for central collisions with the tip (lower barriers) and 
with the side (higher barriers) of the deformed nucleus. From 
Ref. [76]. 



4.2.2 Fusion barriers with a deformed nucleus 

We now consider collisions of a spherical nucleus on a de- 
formed one. In such a case, the barrier depends on the 
orientation of the deformed nucleus at the touching point, 
leading to a wider barrier distribution than the single bar- 
rier case [187,215]. 

Fig. 4.8 shows two examples of experimental barrier 
distributions involving a prolatly deformed heavy target 
[213,214]. Such barrier distributions are usually well re- 
produced in the framework of coupled channel calcula- 
tions [187]. The standard interpretation is that the low 
(resp. high) energy part of the barrier distribution cor- 
responds to collisions with the tip (side) of the deformed 
nucleus. This interpretation has been confirmed with stan- 
dard TDHF calculations [71,76] and density-constrained 
TDHF (DC-TDHF) calculations [70,72,75]. The DC-TDHF 
technique allows for the computation of an energy-dependent 
nucleus-nucleus potential from TDHF trajectories (see Rcfs. 
[216,217,70] for details). The fusion thresholds for the 
two different orientations are represented with arrows in 
Fig. 4.8. This confirms that collisions with the tip of the 
deformed nucleus have lower fusion barriers than collisions 
with the side. We conclude that, in addition to a good re- 
production of the centroids, modern TDHF calculations 
also reproduce the widths of the barrier distributions gen- 
erated by static deformations of heavy targets without any 
adjustment of parameters. 

The case of a light deformed projectile on a heavy 
spherical target has been investigated in Ref. [69] both 
within the TDHF approach and with a coupled channel 
framework [188]. For such systems, the barrier distribution 
is affected by the reorientation of the deformed nucleus in 
the Coulomb field of the target. This induces an anisotropy 
of the orientation axis distribution and results into a fusion 
hindrance at low energies [69,71,75]. Possible experimen- 
tal evidences of this effect have been reported [218]. Note 
that the reorientation is proportional to A sp hericai / A-totai 
and could be neglected in the systems studied in Fig. 4.8. 

Finally, it is worth mentioning that, although TDHF 
calculations help to understand the structure of fusion 
barrier distributions, they would not be able to repro- 
duce their detailed structure due to a lack of quantum 
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Fig. 4.9. Evolution of the expectation value of the dipole mo- 
ment Qd and its conjugated moment Pd in the case of the 
N/Z asymmetric reaction 40 Ca+ 100 Mo at a center-of-mass en- 
ergy of 0.83 MeV/nucleon. 



effects associated to the collective degrees of freedom. In 
particular, the state of the deformed nucleus should be 
a coherent superposition of different orientations in the 
laboratory frame. This quantum coherence is lost at the 
mean-field level. A possible improvement would be to re- 
store this coherence using a time-dependent generator- 
coordinate method (TDGCM) [108], using the orientation 
of the nucleus as a collective coordinate. 



4.2.3 Charge equilibration in fusion 

The charge equilibration process in fusion reactions has 
been investigated with microscopic models such as semi- 
classical [159,219-221] and TDHF [222,68,74,75,89] ap- 
proaches. Charge equilibration occurs when two nuclei 
with different N/Z collide. This difference induces a net 
dipole moment at contact which can oscillate. This isovec- 
tor dipole oscillation is also called preequilibrium GDR, 
as it is a collective motion occuring in the preequilibrium 
stage of the compound system, i.e., before a complete equi- 
libration of its degrees of freedom is reached. 

An example of such preequilibrium dipole motion is 
shown in Fig. 4.9, where the time evolution of the dipole 
moment Qd (proportional to the distance between the 
proton and neutron centers of mass) and its conjugated 
moment Pd (proportional to their relative velocity) are 
shown in the case of the N/Z asymmetric reaction 40 Ca 
+ 100 Mo [74]. Pd and Qd oscillate in phase quadrature. 
They exhibit a spiral in the plot of Pd as a function of 
Q d due to the damping of the dipole vibration. 

It is possible to compute the spectrum of 7 emitted 
by the preequilibrium GDR using laws of classical elec- 
trodynamics. The preequilibrium GDR 7-ray spectrum is 
computed from the Fourier transform of the acceleration 
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Fig. 4.10. (top) 7-spectra measured in the 40 Ca+ 100 Mo 
and 36 S+ 104 Pd reactions at a center-of-mass energy of 0.83 
MeV/nucleon [223]. (middle) Preequilibrium GDR 7-decay 
spectrum obtained from the difference between the two 7- 
spectra in the top. (bottom) Theoretical 7 spectrum computed 
from the preequilibrium dipole moment evolution in Fig. 4.9 
(solid line, purple area). The dotted line represents the first 
chance statistical 7-ray decay spectrum (blue area). Adapted 
from [74]. 
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Such a 7-spectrum is shown in the bottom of Fig. 4.10 
(solid- line, purple area). A comparison with the first chance 
statistical GDR 7-ray decay spectrum is also shown in the 
bottom of Fig. 4.10 (dotted line, blue area). (See Ref. [74] 
for details on the calculations of the latter). We observe 
that the preequilibrium GDR 7 are emitted at a lower en- 
ergy than the statistical component. This is interpreted in 
terms of a large deformation of the compound nucleus in 
its preequilibrium phase [74]. Indeed, the preequilibrium 
dipole motion occurs along the prolate deformation axis 
of the compound system associated to a lower energy (see 
also Fig. 3.4). 

Experimental 7-spcctra are shown in the upper panel 
of Fig. 4.10 for the 40 Ca+ 100 Mo N/Z asymmetric reaction 
and for the 36 S+ 104 Pd reaction which is quasi-symmetric 
in N/Z. Only the first reaction is expected to exhibit a 
preequilibrium dipole motion. Indeed, more 7 are observed 
in this reaction. The difference (purple area in Fig. 4.10- 
top) is interpreted in terms of 7-decay from the preequi- 
librium GDR [223]. Subtracting the two 7-spectra, one ob- 
tains the preequilibrium GDR 7-spcctrum which is shown 
in the middle panel of Fig. 4.10. The energy of the peak 
is in good agreement with the spectrum computed from 
the TDHF response. As mentioned above, this lowering 
of the preequilibrium GDR energy, by comparison to the 
hot GDR, is interpreted as a signature of a strong defor- 
mation of the preequilibrium nucleus. In particular, this 
means that the equilibration of the shape is much slower 
than the charge equilibration. We will discuss charge equi- 
libration times more quantitatively in section 4.5. 

To conclude, we see that the preequilibrium GDR con- 
tains informations on the structure of the preequilibrium 
compound system and, then, on the path to fusion. Here, 
the example of deformation has been discussed. In Rcfs. [68, 
74], the preequlibrium GDR is also used to investigate 
other properties of the path to fusion, such as the role 
of rotation, the coupling with collective shape vibrations, 
the role of mass asymmetry... It is also proposed that the 
decay of a preequilibrium GDR could serve as a cooling 
mechanism increasing the survival probability of the heav- 
iest compound nuclei [221,74]. This effect needs further 
theoretical and experimental investigations. 



4.3 Transfer reactions 



We discussed in section 4.2.1 the interplay between fusion 
and transfer reactions in the 16 O+ 208 Pb system around 
the barrier. In particular, we observed just below the bar- 
rier (see top of Fig. 4.5) an average transfer of two protons. 
We now discuss such transfer reactions in more details. 
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Fig. 4.11. TDHF calculations of central sub-barrier collision 
of two 16 nuclei. The nuclei approach each other, and re- 
separate back-to-back due to the Coulomb repulsion. The evo- 
lution of ap 3 / 2 single-particle wave-function belonging initially 
to the nucleus in the right is shown. After the collision, part of 
this wave function has been transferred to its collision partner. 
Since this reaction is symmetric, a similar transfer occurs from 
the left to the right, and both fragments have the same particle 
number distributions. 
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Fig. 4.12. Transfer in 16 O+ 208 Pb sub-barrier central colli- 
sions. The average number of protons and neutrons of the small 
fragment in the exit channel are plotted as a function of E c . m .. 
The arrow indicates the TDHF fusion barrier. 



Several recent attempts of describing nuclcon trans- 
fer in heavy-ion collisions within the TDHF framework 
have been made in medium mass systems [98,76,42,3,99, 
43] . The TDHF equation describes the evolution of single- 
particle wave-functions. The latter, initially localised within 
one collision partner, may be partially transferred to the 
other fragment during the collision, as illustrated in Fig. 4.11 

In the case of an asymmetric collision, a change of the 
average particle number in the fragments in the exit chan- 
nel is a clear signature that a transfer mechanism occurred 




5 10 

number of nucleons 



Fig. 4.13. Neutron (circles) and proton (squares) number 
probability distributions of the lightest fragment in exit chan- 
nel of a head-on 16 O+ 208 Pb collision at E c . m . = 74.44 MeV 
(top) and 65 MeV (bottom). Adapted from Ref. [99]. 



in the reaction. Figure 4.12 gives the evolution of the ex- 
pectation value of Z and N of the small fragment in the 
exit channel of 16 O+ 208 Pb sub-barrier central collisions. 
At the barrier, ~ 2 protons and no neutron, in average, 
are transferred (the corresponding evolution of the density 
is shown in the top of Fig. 4.5), while at ~ 10% below the 
barrier, Z ~ N ~ 8 is obtained in average, indicating a 
dominance of (in)elastic scattering. We see that the prob- 
ability for proton stripping (transfer from the light to the 
heavy nucleus) is higher than for proton pickup (transfer 
from the heavy to the light nucleus) , while neutron pickup 
is more probable than neutron stripping. This qualitative 
observation is in agreement with experimental data [205, 
100]. 

To get a deeper insight into this transfer mechanism, 
the transfer probabilities are extracted at the TDHF level 
[21,99] thanks to a projection onto a good particle num- 
ber technique 16 applied on the outgoing fragments 17 . It 



16 This technique is standard in beyond-mean-field models for 
nuclear structure when the number of particles is only given in 
average [142]. 

17 One could question this approach as the TDHF wave func- 
tions are not used for the calculation of expectation values of 
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is possible to extract the component of the wave function 
associated to a specific transfer channel using a particle 
number projector onto N protons or neutrons in the x > 
region where one fragment is located at the final time, the 
other one being in the x < region. Such a projector is 
written [225] 
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counts the number of particles in the x > region (0(x) = 
1 if x > and elsewhere). Isospin is omitted to simplify 
the notation. 

The projector defined in Eq. (4.5) can be used to com- 
pute the probability to find N nucleons in x > in the 
final state \d>), 
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where \4>r{6)) = e 1 R \<p) represents a rotation of \<j>) by 
a gauge angle 9 in the gauge space associated to the par- 
ticle number degree of freedom. Note that \4>r(0)) is an 
independent particle state. The last term in Eq. (4.7) is 
then the determinant of the matrix of the occupied single 
particle state overlaps: 



R (0)) =dct(F) 



with 



F tJ =J2 /dr ^f(r)^(r)e^. 



(4.8) 



(4.9) 



The integral in Eq. (4.7) is discretised using 6 n = 27m/M 
with the integer n = 1 • • • M. Choosing M = 300 ensures 
numerical convergence for the 16 O+ 208 Pb system. Fig 4.13 
shows the resulting transfer probabilities at (top) and well 
below the barrier at E c . m . = 65 MeV (bottom). As ex- 
pected from the average values (see Fig. 4.12), the most 
probable channels arc Z = 6 and TV = 8 at the barrier, 
and Z = N = 8 with a small probability of neutron pickup 
or proton stripping (of the order of 10~ 2 ) well below the 
barrier. 

A standard representation of experimental sub-barrier 
energy transfer data is to plot transfer probabilities as a 
function of the distance of closest approach Rmin between 
the collision partners [226]. R m in is computed assuming a 
Rutherford trajectory [227]: 



R mm = ZiZ 2 e 2 [l + cosec(0 c . m ./2)]/2£ c . r 



(4.10) 



where 



is the center of mass scattering angle. 



one-body operators. In particular, the width of the distribu- 
tions should be underestimated [23]. However, as we will see 
in section 4.4, the TDHF and BV widths are similar for non 
violent collisions such as sub-barrier transfer, justifying a pos- 
teriori this approach. 
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Fig. 4.14. Proton number probability as function of the dis- 
tance of closest approach in the small outgoing fragment of 
the 16 O+ 208 Pb reaction. TDHF results are shown with lines. 
Experimental data (open symbols) are taken from Ref. [100]. 



A comparison of TDHF calculations with recent data 
from Ref. [100] is shown in Fig. 4.14 for sub-barrier one 
and two-proton transfer channels in the 16 O+ 208 Pb re- 
action. We see that TDHF overestimates the one-proton 
transfer probabilities and underestimates the 2-proton trans- 
fer channel. This discrepancy is interpreted as an effect 
of pairing interactions [99,100]. Indeed, paired nucleons 
are expected to contribute to two-nucleon transfer chan- 
nels. As a result, the two-nucleon (resp. single-nucleon) 
transfer probability increases (decreases). Note, however, 
that for R m i n > 13 fm, the TDHF calculations reproduce 
reasonably well the sum of one and two-proton transfer 
channels 18 [99]. 

These studies emphasize the role of pairing interac- 
tions in heavy-ion collisions. The recent inclusion of pair- 
ing interactions in 3-dimensional microscopic codes [17, 
114,65] gives hope in our ability to describe such data 
with more details in a near future. Finally, microscopic 
theories should also be used to determine, not only the 
transfer probabilities, but also the excitation energy of the 
final fragments. In particular, this would be of great help 
to understand how transfer reactions induce an energy 
dissipation, possibly hindering fusion at deep sub-barrier 
energies [105,100]. 



4.4 Deep-Inelastic Collisions 

Deep- Inelastic Collisions (DIC) occur essentially well above 
the barrier. The main characteristics of DIC exit channels 



A strong damping of the initial kinetic energy, 
Large fluctuations of the fragment proton and neutron 
numbers around their initial value, 
An angular distribution of the fragments following a 
1 / sin B c . m . behaviour. 



18 For R m i„ < 13 fm, sub-barrier fusion, not included in the 
TDHF calculations, reduces transfer probabilities [99,100]. 
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Fig. 4.15. Density evolution for a 40 Ca+ 40 Ca collision at 
E c . m . = 128 MeV and L — 60ft. Each snapshot is separated 
by 1.5 zs. 



The last point is due to a large orbiting of the fragments at 
contact. All 6 c . m . are then equiprobable in DIC. However 
the emission is not isotropic as it occurs essentially in the 
collision plane. As a result, the differential cross-sections 
for DIC events obey 

da da 

— = z7T sin u c . m . — — ~ constant. 

du crn die 

This leads to the 1/ sin# cm . behaviour, as 
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Early TDHF calculations were able to reproduce frag- 
ment kinetic energies, mean masses, and scattering angles, 
but they failed to reproduce the observed large fluctua- 
tions of the fragment Z and N distributions [21,22]. This 
failure was seen as a necessity to include a collision term in 
the dynamics. However, Balian and Veneroni showed that 
the TDHF theory was in fact not optimised for the predic- 
tion of such fluctuations [20]. Instead, one should use the 
prescription given in Eq. (2.40) which is derived from their 
variational principle and optimised for one-body fluctua- 
tions [24]. 

The BV prescription was used to compute particle 
number fluctuations in DIC in Refs. [28,27,31,32,228], 
and within a semi-classical approximation in Ref. [29]. 
Comparing with standard TDHF calculations, these works 
showed, indeed, an increase of the fluctuations. However, 
the first realistic calculations allowing for a direct com- 
parison with experimental data 19 was only recently per- 
formed in Ref. [32] . We now discussed the main results of 
this paper. 

40 Ca+ 40 Ca collisions at E c . m . = 128 MeV (~ 2.4 times 
the barrier height) have been studied with the tdhf3d 
code. Figs. 4.15, 4.16, and 4.17 show density plots obtained 



19 The authors of Ref. [27] claim that their calculations are in 
good agreement with experimental data. However, they only 
compute the fluctuations of A — N + Z and compare with 
measured fluctuations of Z. In addition, their calculations do 
not include the spin-orbit interaction. They are performed at 
angular momenta leading to fusion when spin-orbit terms are 
included (see discussion in Ref. [31]). 



Fig. 4.16. Same as Fig. 4.15 with L = 70ft. 



Fig. 4.17. Same as Fig. 4.15 with L = 80ft. 





Fig. 4.18. Trajectories of the centers of mass of the fragments 
in 40 Ca+ 40 Ca collisions at E c . m . = 128 MeV. 



with angular momenta L=60, 70, and 80 (in units of ft), re- 
spectively Three reaction mechanisms are observed: cap- 
ture (at L = 60), DIC (at L = 70), and a partially damped 
collision (at L = 80). In particular, the L = 70 case leads 
to an orbiting trajectory (see middle panel of Fig. 4.18) 
which is characteristic of a DIC. This trajectory also cor- 
responds to a strongly damped collision, as can be seen 
from the upper panel of Fig. 4.19 where the total kinetic 
energy loss (TKEL) is plotted as a function of the initial 
angular momentum. Indeed, around L ~ 70, the TKEL is 
~ 60 — 70 MeV. These values are slightly below the Viola 
systematics [229] which predicts TKELyioia — 76 MeV, 
indicating that these collisions are almost fully damped. 

The lower panel of Fig. 4.19 shows the evolution of the 
fluctuations az of the number of protons 20 in the outgo- 
ing fragments as a function of L. As expected 21 , the BV 
prescription [Eq. (2.40)] leads to larger fluctuations than 
standard TDHF fluctuations obtained from Eq. (2.37). 
This is particularly true for DIC. However, TDHF and 
BV prescriptions converge at large L, e.g., for L > 90 
where the TKEL is less than 5 MeV, indicating a domi- 



20 Neutron fluctuations are very close to the proton ones for 
this system. This is due to the fact that the collision partners 
are N = Z nuclei. 

21 See early applications of the BV prescription for fluctua- 
tions where greater fluctuations than their TDHF counterparts 
were obtained [27,28,31]. 
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Fig. 4.19. Properties of the exit channel of the 40 Ca+ 40 Ca 
collisions at E c . m . = 128 MeV as a function of the initial angu- 
lar momentum L. (top) Total kinetic energy loss from TDHF. 
(bottom) TDHF (dashed line) and BV (circles) fluctuations of 
Z, and BV correlations between N and Z (triangles) of the 
outgoing fragments. 



nance of quasi-elastic scattering. We conclude that TDHF 
calculations of particle number distributions give reason- 
able estimates for such non- violent collisions. In particu- 
lar, this justifies the calculation of transfer probabilities 
for sub-barrier collisions with TDHF (see section 4.3). 

Fig. 4.19 also shows an application of the BV prescrip- 
tion for the correlations between N and Z distributions. 
These correlations are determined from Eq. (2.38). We ob- 
serve increasing correlations with decreasing L. In partic- 
ular, in DIC, the BV prescription predicts 22 <jn — <jz — 
onz-i while, for quasi-elastic scattering, correlations are 
negligible. 

Let us recall the signification of such correlations. Un- 
corrected distributions mean that the probability to trans- 
fer z protons is independent of the probability to transfer 
n neutrons, i.e., P(z,n) = P(z)P(n). On the other side, 
strongly correlated distributions mean that if we measure 
n (rcsp. z), then we know what is z (n). This is the case, 
for instance, if all fragments have N = Z. In this case, we 
would have P(z,n) = P(z)S n , z = P(n)S n , z - The reality is 
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Fig. 4.20. Comparison between BV (solid line) and TDHF 
(dashed line) predictions of azz for damped events (see text) 
as a function of 6 c . m . with data (circles) from [230]. 



usually in between and, assuming Gaussian distributions, 
we have 
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where 



= — iLZ -. The case p = means no correlations 
between N and Z distributions, while the limit \p\ 



1 



(tat is not shown in Fig. 4.19 for the clarity of the figure. 



corresponds maximum correlations. 

In the calculations of 40 Ca+ 40 Ca at E c , m . = 128 MeV 
shown in Fig. 4.19, we have unz — &n — <Jz for DIC 
(L < 80). This means that N and Z distributions of the 
fragments are strongly correlated in DIC. However, quasi- 
elastic reactions (L > 90) have <tn — oz S> <jnz, meaning 
almost independent N and Z distributions in this case. 

We now compare these results with the experimental 
data of Roynette et al. [230]. The angle c . m . between the 
fragments in the outgoing channel and the collision axis 
have been determined for each L. Fig. 4.20 shows theo- 
retical and experimental evolutions of the charge fluctu- 
ations for damped events (defined, as in Ref. [230], by 
a TKEL > 30 MeV) as a function of 6 c . m .. Due to or- 
biting, only DIC are expected to emit fragments at large 
angles, and the experimental plateau at c . m . > 50 deg is 
then attributed to DIC. We see that TDHF fluctuations 
underestimate experimental results at all angles, except 
at very forward angles where quasi-clastic reactions dom- 
inate. The results of the BV prescription are in better 
agreement, although they still underestimate the experi- 
mental data. This is probably due to fusion-fission events 
(not included in the calculations) leading to large fluctu- 
ations and, to a less extent, to the cooling down of the 
fragments by nucleon emission [32] . 

In addition, calculations with exotic nuclei have been 
performed to study the role of isospin asymmetry in the 
entrance channel on the Z and N fluctuations and corre- 
lations [32]. The 80 > 92 Kr+ 90 Zr systems have been investi- 
gated at a beam energy E/A = 8.5 MeV. Charge equi- 
libration is observed in the 92 Kr+ 90 Zr due to an initial 
N/Z asymmetry, inducing larger correlations between TV 
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and Z distributions, while fluctuations are only slightly 
affected. This increase of the correlations between N and 
Z distributions due to charge equilibration could be tested 
in DIC experiments with exotic beams [231]. Note, finally, 
that other fluctuations should also be computed with the 
BV approach, such as the width of TKE distributions of 
the fragments. 



4.5 Fusion and quasi-fission in heavy systems 

In section 4.2, we restricted the study to Z1Z2 < 1600 sys- 
tems for which no fusion hindrance is usually observed 23 . 
We now investigate the reaction mechanism in heavy sys- 
tems with possible fusion hindrance. We first illustrate the 
fusion hindrance with TDHF calculations of fusion thresh- 
olds [157,202]. Then, we present recent results on a study 
of the quasi-fission process with the TDHF approach. 



4.5.1 TDHF calculations of fusion hindrance 

As an example, we consider the 90 Zr+ 124 Sn system which 
has a charge product Z\Z- 2 = 2000, and, then, is expected 
to exhibit a fusion hindrance. The proximity model [235] 
predicts a barrier for this system V prox ' ~ 215 MeV. 
Fig. 4.21 shows the TDHF evolution of the relative dis- 
tance between the fragments as a function of time for 
central collisions at different energies [157]. We see that 
the system encounters a fast re-separation at the energy 
of the barrier predicted by the proximity model. Long 
contact times possibly leading to fusion are observed at 
E c . m . > 240 MeV. 

Density profiles for this system at E c . m . = 235 MeV 
are shown in Fig. 4.22. A rapid neck formation is observed. 
However the system keeps the shape of two fragments 
in contact during ~ 1400 fm/c before re-separation in 
two fission-likc fragments. This reaction mechanism differs 
from fusion followed by statistical fission as a compound 
system is not formed in the present case. In particular, 
the system is expected to keep the memory of its entrance 
channel. Such a process is called quasi-fission. 

Similar calculations have been performed for other sys- 
tems with ZxZi > 1600 [202]. The TDHF fusion thresh- 
olds arc shown in Fig. 4.23. A comparison with the in- 
teraction barriers predicted by the proximity model [235] 
shows that dynamical effects included in TDHF induce 
a strong increase of the fusion threshold, in particular for 
the more heavy and symmetric systems. The order of mag- 
nitude of the additional energy needed to fuse is similar to 
the one predicted with the extra-push phenomenological 
model [232]. 



This threshold is empirical. Based on his extra-push model, 
Swiatecki proposed an effective fissility, depending on both 
charges and masses of the nuclei, above which extra-push en- 
ergy is needed to fuse [232]. Note that this should not be con- 
founded with a threshold for quasi-fission which, in fact, may 
occur in lighter systems. Indeed, quasi-fission has been ob- 
served in, e.g., 16 0, 32 S+ 238 U [214,233,234], and even in lighter 
systems such as 32 S+ 208 Pb [234]. 
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Fig. 4.21. Distance between the centers of mass of the frag- 
ments as a function of time in head-on 90 Zr-|- 124 Sn collisions 
for different center of mass energies. From Ref. [157]. 




Fig. 4.22. Density profile in the 90 Zr+ 124 Sn head-on collision 
at Ecm. = 235 MeV. From Ref. [157]. 
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Fig. 4.23. TDHF fusion thresholds for several heavy systems 
are compared with the proximity barrier [235] and with results 
from the extra-push model [232]. 



These results motivate a deeper study of the origin of 
the fusion hindrance in such heavy systems. In particu- 
lar, the quasi-fission mechanism which is mostly responsi- 
ble for this hindrance is studied in section 4.5.2 with the 
TDHF approach. 



4.5.2 TDHF calculations of the 40 Ca+ 238 U reaction 

The previous studies showed the importance of the quasi- 
fission process as a mechanism in competition with fusion, 



Cedric Simenel: Nuclear Quantum Many-Body Dynamics 



33 



Fig. 4.24. Initial condition of a 40 Ca+ 238 U collision visualised 
with the SDVision code [236]. Two isodensities are shown (yel- 
low grid and blue area). A projection of the density is also 
shown on the z = plane. 
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Fig. 4.25. Snapshots of the TDHF isodensity at half the sat- 
uration density in the 40 Ca+ 238 U system for different initial 
orientations and E c . m .. 



hindering the formation of heavy systems. A quantum 
and microscopic theoretical framework able to describe 
properly the quasi-fission properties would be of great im- 
portance to get a deep insight into the interplay between 
structure properties and this mechanism. 

Here, our goals are to show that quasi-fission may ap- 
pear in the outgoing channel within the TDHF approx- 
imation, and that TDHF calculations can then be used 
to predict quasi-fission properties. An example of initial 
condition of the 40 Ca+ 238 U collision is shown in Fig. 4.24 
where we can clearly see internal structures in the 238 U 
nucleus. 
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Fig. 4.26. TDHF calculations of the mass of the fragments 
(top) and of the quasi-fission time (bottom) in 40 Ca+ 238 U cen- 
tral collisions as a function of the center of mass energy (di- 
vided by the proximity barrier [235]). For quasi-fission times 
larger than 23 zs, only a lower limit is given. Two differ- 
ent orientations of the 238 U are considered (see inset). From 
Ref. [237]. 



Examples of density evolutions obtained with the tdhf3d 
code are shown in Fig. 4.25 for the 40 Ca+ 238 U system 
around the barrier [237] . These configurations lead to quasi- 
fission as we observe an important multi-nuclcon transfer 
from the heavy fragment toward the light one. We also 
see that the mass equilibration (i.e., the formation of two 
fragments with symmetric masses) is not complete and 
may depend on the initial conditions. 

Fig. 4.26 presents results for 40 Ca+ 238 Ca central colli- 
sions. Final fragment masses (top) and quasi-fission times 
(bottom) are shown for two different orientations of the 
238 U. We see that all the calculations with the 238 U defor- 
mation axis aligned with the collision axis lead to a quasi- 
fission with partial mass equilibration and quasi-fission 
times smaller than 10 zs. Shell effects may affect the final 
outcome of the reaction by favouring the production of 
fragments in the 208 Pb region. In particular, this orienta- 
tion never leads to fusion, while the other orientation pro- 
duces long contact time above the barrier which may be 
associated to fusion. We also see that longer quasi-fission 
times lead to larger mass equilibrations. Calculations of 
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Fig. 4.27. Nucleus-nucleus potentials of the 238 TJ+ 238 U sys- 
tem obtained with HF ground-states and the frozen approxi- 
mation for different orientations of the actinides as a function 
of the distance between their centers of mass. 
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Fig. 4.28. The adiabatic potential is computed with the ev8 
code [240] as a function of the total quadrupole moment. Cal- 
culations are performed with (HFB) and without (HF) pairing 
residual interaction. 



non-central °Ca+ 238 U collisions are ongoing in order to 
compare with experimental data. 



4.6 Actinide collisions 

The collision of actinides form, during few zs, the heavi- 
est nuclear systems available on Earth. In one hand, such 
systems are interesting to study the stability of the QED 
vacuum under strong electric fields [238,239,102]. In the 
other hand, they might be used to form neutron-rich heavy 
and super-heavy elements via multi-nuclcon transfer reac- 
tions, as discussed in section 4.1.2. 

As we saw in section 4.5, there is a link between the col- 
lision time and the amount of transferred nucleons [197]. 
It is then important to optimize collision times in or- 
der to favour the formation of heavy systems. It was ini- 
tially believed that the potentials between two actinides 
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Fig. 4.29. Snapshots of the isodensity at half the saturation 
density in 238 tj+ 238 TJ central collisions at E c . m . = 900 MeV 
from TDHF calculations. Snapshots are given at times t — 0, 
1.5, 2.7, and 4.2 zs from top to bottom. From Ref. [102]. 



have a barrier (and then a pocket) [241], leading to pos- 
sible long contact times at energies close to the barrier. 
Such barriers are predicted by frozen models, as shown 
in Fig. 4.27 where the frozen potentials have been com- 
puted with two 238 U HF ground states with different ori- 
entations. However, Berger and collaborators showed with 
constrained HFB calculations of the composite system 
of 238 U_|_238tj that there is, in fact, no barrier in their 
(adiabatic) nucleus- nucleus potential [242]. This result is 
confirmed with modern Skyrme-HF calculations using the 
EV8 code [240], as shown in Fig. 4.28. 

The dynamics of actinide collisions have been studied 
with the TDHF approach [243, 102, 103, 228] and the quan- 
tum molecular dynamics (QMD) formalism [244,245]. These 
microscopic studies complement other works with a quan- 
tum master equation [246], the dinuclear system (DNS) 
model [247, 248] , multidimensional Langevin equations [249] , 
and the constrained molecular dynamics model [250]. 

Fig. 4.29 shows zeptosecond snapshots of the density 
obtained in TDHF calculations of 2 38tj-|_238tj central col- 
lisions at -E c . rn . = 900 MeV. We see that the initial ori- 
entation of the nuclei plays a crucial role on the reaction 
mechanism [102]. For instance, wc observe the formation 
of a third fragment in the left column. A net transfer is 
also obtained in the middle column. Indeed, integration of 
proton and neutron densities in each reactant indicates an 
average transfer of ~ 6 protons and ^11 neutrons from 
the right to the left nucleus. In this case, transfer occurs 
from the tip of the aligned nucleus to the side of the other. 
This configuration is then expected to favor the formation 
of nuclei heavier than 238 U. 

Pursuing this idea, a new "inverse quasi- fission" mech- 
anism has been identified [103]. The inverse process of 
quasi-fission corresponds to a transfer of nucleons from 
the light collision partner to the heavy one. Such a mech- 
anism may occur in actinide collisions due to shell effects 
in the 208 Pb region [251,249]. The mechanism proposed 
in Ref. [103] is different (but may be complementary). It 
occurs for specific orientations of the actinides, where the 
tip of the lighter one is in contact with the side of the 
heavier one. 
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An example of inverse quasi-fission reaction is shown in 
the right panels of Fig. 4.30, where a 232 Th+ 250 Cf collision 
produces a 265 Lr fragment at the end of the calculation 24 . 
The 265 Lr heavy fragment indicates the average N and Z 
of a distribution. The fluctuations and correlations of these 
distributions have been computed with the BV prescrip- 
tions [Eqs. (2.40) and (2.38)] in Rcf. [228]. Fig. 4.30(left) 
shows the resulting probabilities assuming Gaussian dis- 
tributions of the form given in Eq. (4.11). We see that 
many j3— stable and neutron-rich transfermium primary 
fragments could be produced thanks to this inverse quasi- 
fission mechanism. In particular, these nuclei are more 
neutron-rich than those formed in fusion-evaporation re- 
actions. This inverse quasi-fission process needs further 
studies to determine the role of the shape, orientation, 
shell closures, mass and charge asymmetries, beam energy, 
and angular momentum on the transfermium production 
yields. Associated cross-sections need to be determined 
to estimate the experimental possibility of neutron-rich 
transfermium and SHE productions. 

Only few early experimental data exist on heavy ele- 
ment production in actinide collisions, with chemical iden- 
tification techniques leading to limited sensitivity in terms 
of cross-section measurements. For example, the smallest 
cross-sections obtained for the production of transfermi- 
ums (in this case, mendelevium isotopes with Z = 101) 
were of the order of ~ 20 nb in the 238 U+ 248 Cm reac- 
tion [252], while modern fusion experiments using spec- 
trometers for the fragment separation have measured cross- 
sections down to 30 fb [195], i.e. ~ 6 orders of magni- 
tude smaller. The use of modern experimental techniques 
and equipments to identify multi-nucleon transfer prod- 
ucts with actinide collisions is then expected to be very 
fruitful in terms of transfermium studies [253] . 

A non exhaustive list of relevant heavy nuclei to be 
searched for is given below: 

— SHE in the island of stability (e.g., 298 114). 

— 2 6i-263p m w hi cn have a predicted /3-decay. If observed, 
they would be the heaviest nuclei with this decay mode. 

— 264 Fm could undergo symmetric spontaneous fission to 
form two doubly magic 132 Sn. 

— From 271 Rfio4 to 284 Mti g, /3— stable nuclei and their 
neighbours are predicted to decay by spontaneous fis- 
sion. Verifying this prediction would constrain the mod- 
els predicting a spontaneous fission region in the decay 
path of r-process progenitors with mass ~ 300 [189]. 
The existence of the latter is crucial to determine if 
long-lived SHE can be produced in the r-process. 

— 264 Noio2 has 162 neutrons which corresponds to a de- 
formed shell gap in the Z = 108 region. The question 
of the robustness of the N=162 shell gap away from 
Z=108 could be answered by measuring the mass and 
life-time of 264 Noio2 and its neighbours. 

24 This fragment is not exactly a primary fragment as at the 
final time of the calculation, about 3 neutrons have been emit- 
ted in the entire system (see Fig. 4 of [103]). However, it is 
still excited (its excitation energy may be estimated from the 
TKE) and may cool down by the emission of 1-2 additional 
neutrons or, of course, by fission. 



— 266 Rfio4 and its neighbours should be studied to "fill 
in" the blank spot between cold and hot fusion decay 
chains. 

— 290 114, 287 113, 286 112 should be produced to confirm 
the increase of life time with neutron number attributed 
to the closeness of the island of stability 25 . 

— Elements Z = 109— 111 still need chemical character- 
isation [192]. Present chemical techniques tag the iso- 
topes thanks to their a— decay and Ti/ 2 > 0.5 s is re- 
quired. 278 Mtio9 and 282 Rgm are good candidates. Ds 
(Z — 110) element is more problematic as no isotope 
decaying by a— emission with T\i% > 0.5 s is known. 
Decay modes of 280 Dsuo should then be measured. If 
the latter does not match experimental requirements 
for chemistry studies, isotopes on the other side of the 
spontaneous fission region (see Rcf. [189]), i.e. 287 Ds 
and its neighbours, might be considered as they arc 
expected to decay by a or /3~ emission. 

— 2 115, which has T\/2 — 0.7 s is the next milestone in 
chemical studies of SHE. 



4.7 Conclusions and perspectives 

This section was devoted to the study of heavy-ion reac- 
tion mechanisms in order to investigate quantum dynami- 
cal effects in complex systems and to predict the formation 
of new nuclei. 

The TDHF theory was shown to reproduce energy 
thresholds for fusion for systems spanning the entire nu- 
clear chart. The effect of deformation and transfer on fu- 
sion barriers is well treated. The competition with quasi- 
fission leading to a fusion hindrance in heavy systems is 
also included. The possibility to study the quasi-fission 
mechanism with a fully microscopic quantum approach is 
promising. It will help to understand the strong fusion 
hindrance in quasi-symmetric heavy systems. 

It was also shown that the charge equilibration process 
in N/Z asymmetric collision affects the reaction mech- 
anism: excitation of a preequilibrium GDR in the com- 
pound system (which can be used to study the path to fu- 
sion), enhancement of proton stripping and neutron pickup, 
and correlations between fragment N and Z distributions 
in DIC (which should be sensitive to the symmetry en- 
ergy). 

Actinide collisions have been investigated both within 
the TDHF approach and with the BV prescription. A new 
inverse quasi-fission mechanism associated to specific ori- 
entations was found in numerical simulations. This mech- 
anism might produce /3— stable and neutron-rich heavy 
and super-heavy nuclei. A systematic investigation of this 
effect (cross-sections, angular and energy distributions) 
is mandatory to help the design of future experimental 
equipments dedicated to the study of fragments produced 
in actinide collisions. 



25 277 112 and 285 112 have T 1/2 = 0.7 ms and 34 s measured 
half-lives, respectively. A gain of 8 neutrons increases the life- 
time by 5 orders of magnitude! 
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Fig. 4.30. (right) Snapshots of the isodensity at half the saturation density in 232r ph+ 250 Cf central collisions at E c . m . = 
916 MeV. (left) Gaussian distributions of N and Z of the heavy outgoing fragment with widths and correlations computed with 
the BV prescription (linear color scale). The purple line shows the expected /3-stability line. 



Quantum effects at the single-particle level are well 
treated in the TDHF approach, allowing for realistic pre- 
dictions of independent particle transfer probabilities us- 
ing a particle-number projection technique. A comparison 
with experimental data shows the importance of pairing 
correlations on transfer. The latter could be studied with 
3-dimensional TDHFB codes in the future. In particular, 
one could answer the question on the origin of these pair- 
ing correlations: Are they present in the ground-states or 
are they generated dynamically during the collision? 



A strong limitation of the TDHF approach is that it 
does not allow for quantum tunnelling of the many-body 
wave function. As a result, sub-barrier fusion cannot be 
studied. Recently, suppression of tunnelling was observed 
at deep sub-barrier energies, which could not be repro- 
duced by standard coupled channel calculations [254,105, 
255] . This is one of the biggest puzzles in low-energy heavy 
ion collision physics, attracting the interest of physicists 
from various fields including astrophysics, where changed 
quantum tunnelling rates can have drastic effects [256]. 
In the future, microscopic models should be developed to 
study such low-energy reactions including the tunnelling 
of the many-body wave-function. A possible approach is 
to use an extension of the TDHF theory based on a path 
integral technique and on the stationary phase approxi- 
mation [45] (sec discussion in section 2.7). 



5 Conclusions 

Nuclei are ideal to investigate fundamental aspects of the 
quantum many-body problem. They exhibit collective mo- 
tions built from coherent superpositions of the states of 
their constituents. Examples range from collective vibra- 
tions to the formation of a compound system in collisions. 
These features are common to other composite systems 
(atomic clusters, molecules, Bose-Einstein condensates...). 
Their study in nuclear systems is obviously part of a wider 
physics field. 

The Balian-Vcncroni variational principle offers a pow- 
erful approach to the many-body quantum dynamics. Both 
the observable and the state of the system are considered 
as variational quantities. In the limit of independent par- 
ticles, one obtains the time-dependent Hartree-Fock for- 
malism for the expectation of one-body observables, and 
the time-dependent RPA for their fluctuations and corre- 
lations. 

Different studies of nuclear dynamics, from collective 
vibrations to heavy-ion collisions have been presented in 
this review. A particular attention was devoted to the in- 
terplay between collective motions and internal degrees of 
freedom within a unified theoretical description. 

We showed that strongly interacting systems such as 
nuclei could exhibit collective vibrations in the continuum, 
and that their direct decay could be used to infer their 
microscopic structure. We also questioned the harmonic 
nature of these vibrations. In particular, we identified a 
source of anharmonicity from the coupling between differ- 



Cedric Simenel: Nuclear Quantum Many-Body Dynamics 



37 



ent vibrational modes. Nuclei are known to exhibit super- 
fluidity due to pairing residual interaction and it is natural 
to wonder if the pairing field could also vibrate. We then 
extended the theoretical approach to study such pairing 
vibrations. 

Large amplitude collective motions were investigated 
in the framework of heavy-ion collisions. We described the 
mechanism of the formation of a compound system. In 
particular we discussed how fusion is affected by the in- 
ternal structure of the collision partners. We then inves- 
tigated the other mechanisms in competition with fusion, 
and responsible for the formation of fragments which dif- 
fer from the entrance channel: transfer reactions, deep- 
inelastic collisions, and quasi-fission. We finally studied 
actinide collisions forming, during very short times of few 
zeptoseconds, the heaviest nuclear systems available on 
Earth. In particular, we identified a new reaction mecha- 
nism occurring in these collisions which could be used to 
study the upper part of the nuclear chart. 
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APPENDICES 

A TDHF approach from the standard 
variational principle 

Consider a state of A*" independent particles described by 
the Slater determinant \<p). The action defined in Eq. (2.2) 



reads 

S = S to , tl [4>] 



\ m)\(ihj t -H 



lit 



where H is the Hamiltonian of the system. In such a state, 
every M-body density matrix (1 < M < N) is simply 
expressed as a function of the one-body density matrix 
p (see appendix C of Ref. [3] ) . Then, all the information 
on the system is contained in p. The expectation value of 
the Hamiltonian on the state \<p) may then be written as 
a functional of p: E[p] = (<f>\H\(f)). In addition, we have 



JV 



(A.2) 



and, then, 



N 



s = jf dt UhJ2(vi\j t \vi) - E[ P (t) 



i=l 

dt (i h J2 J dx V*^'*) ^Pi{x,t) - E\p{t)\ 



(A.3) 



where x = (rsq) describes all the single-particle degrees of 
freedom (position r, spin s and isospin q). We introduced 
the notation Jdx = J2 qs fd r - 

The variational principle reads SS = 0. The variation 
must be done on each independent variable. Here, these 
variables are the real part (p^ c and the imaginary part ip 1 ™ 
of each occupied single particle state <p a . We must then 
consider 



SS 



= 



and 



SS 



S^(x,t) 



(A.4) 



for each a € {1...N}, for all t such as to < t < t\, and for 
all x. 

However, the calculation is more straightforward if we 
use (p a and <p* a as independent variables instead of tp^ c and 
<p 1 ™- We can then consider the variations over <p and over 
tp* independently. Note, however, that we loose the prop- 
erty "v?* being complex conjugated of ip" that we should 
restore later. 

Equations (A.4) are then replaced by 



SS 



= 



and 



SS 



8<p a (x,t) 

The variation over ip* gives 

SS fc d . . 
8<p* a {x,t) dt Jt 



8<p* a {x,t) 



= 0. (A.5) 



\, 5E[p{t')] 
o Sp* a (x,t) 



The functional derivative of E can be re-written thanks 
to a change of variable 



SEW)] = f dvdv , SEW')} 5p(y,y';t>) 
Sip* a (x,t) J V V Sp(y,y';t') S<p* a (x,t) 



(A.7) 
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Using 



sp(y,y';t) 

Sip* a (x 7 t) 



tp a (y,t')5(y'-x)S(t-t') (A. 



and noting the singlc-particlc Hartrcc-Fock Hamiltonian 
h with matrix elements 



h(x,y\t) 



SE[p(t)} 



(A.9) 



we get the TDHF equation for the set of occupied states 



ih ^ a {x,t) = jdy h(x,y;t) tp a (y,t) 



(A.10) 



The variation over tp gives, after integrating by part 
the term with the time derivative, 



SS 



5tp a {x,t) 5tp a (x,t) J ta 



dt' 



ih ~& \ Y1 J dy tfsiyJ) <pp(vJ) 

(A.ll 



The first term in the r.h.s. cannot be replaced a priori by 
Tr[p] = N because we considered tp and tp* to be indepen- 
dent variables. Let us note the variation of tp at a specific 
"point" of the Hilbcrt space at time t and affecting only 
the state tp a . This variation can be defined as 



B Fluctuations and correlations of one-body 
observables with the Balian-Veneroni 
variational principle 

In this appendix, we detail the derivation of the Balian- 
Vcncroni prescription for fluctuation and correlation of 
one-body observables, given in Eqs. (2.40) and (2.38), re- 
spectively. We consider the particular case where the in- 
dependent particle state of the system is a pure state de- 
scribed by a unique Slater determinant at all time. The 
more general derivation with mean-field states of the form 
given in Eq. (2.32) can be found in Refs. [25,31]. 

As mentioned in section 2.4, fluctuations and correla- 
tions of the one-body observables 



a/3 



(B.1) 



(B.2) 



at the final time t\, and using Eq. (2.29) for small The 
expectation value of A\ at t\ is obtained from 



are obtained by evaluating 



A\ = cxp — 




(A 1 )(t 1 )=Tr A(h)D(t 



(B.3) 



where the observable is assumed to be a time-dependent 
operator A(t) with the boundary condition A(t\) = A\. 
The state of the system is assumed to be known at the 
initial time t^ 1 with the boundary condition D(to) = Dq. 
As a result, the variations of the observable and the state 
obey 

SA(tx) = SD(t ) = 0. (B.4) 



5 axt ipp{y,t') = Aip f(t) 5(t-t') d(x-y) S af3 . (A.12) 

Using the usual definition of the functional derivative, we 
may write 

= 4™o//*' ih ~h I dy ^ {y,t ' ] ^ t ^ (j/, *' ) ) 

= [/(*) 5(t-t')v* a (x,t%L t0 . (A.13) 

We have to choose /(io) = f(ti) = in oder to cancel 
this term at t = to or t = t\. It is equivalent to forbid 
variations of ip at to and t\. As a result, Eq. (A.ll) leads 
to the complex conjugated of the TDHF equation (A. 10), 
restoring the fact that tp and tp* are complex conjugated. 
This last property is also necessary for energy and particle 
number conservations. 



Variational space and parametrisation of variational quan- 
tities 

The BV variational principle is solved by requiring the sta- 
tionarity of the action-like quantity J defined in Eq. (2.12), 
or, equivalently, in Eq. (2.14). In the present application, 
the variational space for the observable is restricted to ex- 
ponential of one-body operators. The observable is then 
parametrised as 

A(t) = exp[-L(t)} (B.5) 

with 

L(t)=Y,L a[i {t)aia p (B.6) 

a/3 

The mean-field density-matrix is parametrised as 

D(t) = cxp[-m(t) - M{t)\ (B.7) 

with 

M(t)=J2 M Mt)ata - (B.8) 

a/3 
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As discussed in section 2.4, the case of a pure Slater de- 
terminant corresponds to the limit where the eigenvalues 
of the matrix M(t) tend to ±00 and m(t) — > +00 in such 
a way that z(t) = TrD(t) = 1 [39]. 

According to Eqs (C.l) and (C.17), we have the fol- 
lowing relationships 



Pit) 



and 



1 + e M (*) 
m[t) =trln(l + e- M W). 



(B.9) 
(B.10) 



where p is the one-body density matrix with elements 

p a p(t) = Tr[b(t)tfph a \. (B.ll) 
The inverse relationships read 

m(t) = -tr[ln(l - p(t))] and (B.12) 



-M(t) 



-tr[ln(l - p(t))] and 
P(t) 



1 - p(t) 



(B.13) 



As expected, we see that the Slater determinant can be 
entirely parametrised by p(t). 

By analogy, we introduce a similar parametrisation of 
the observable using Eqs. (C.l) and (C.12). We get 



a(t) 

y(t) 



1 



1 + e L W 



and 



exp 



trln(l 



with the inverse relationships 

lny(t) = -tr[ln(l - a(t))} and 

-£(*) = °(t) 

l-tr(t)" 



(B.14) 
(B.15) 

(B.16) 
(B.17) 



Eqs. (2.12) and (2.14) contain the products AD and 
DA. It is convenient to define a similar parametrisation 
of these quantities, using the fact that the product of the 
exponential of one-body operators is also an exponential 
of a one-body operator 



A(t)D{t) = e- L ^e- M W 

b(t)A{t) = e -m) e -L{t) 



~ L ' W and (B.18) 



(B.19) 



-M'(t) 



L' and M' are one-body operators which are parametrised 
as 



with 



L'{t)=m{t)+Y J L' a ^t)al l ap and 

a/3 

M'{t)=m{t) + Y J M' a0 {t)aiap, 

a/3 



e -L'(t) = e -L{t) e -M{t) and 
e -M'(t) = e -M(t) e -L(t)_ 



(B.20) 
(B.21) 



(B.22) 
(B.23) 



By analogy with Eq. (B.ll), we introduce two new matri- 
ces, p' and a' ', with elements 



Tv[D{t)A{t)a} p a 
TT:[A{t)b{t)a\d 



and 



and obeying the following relationships: 
1 1 

P'(t) = 



a'(t) = 



l + e M'(t) l + e L(t) e M(t) 
1 1 



(B.24) 
(B.25) 

and (B.26) 
(B.27) 



l + e £'(t) i + e M(t) e L(t)- 
Using Eqs. (B.19) and (C.12), the normalisation factor 

u{t) = Tr[D(t)A(t)] (B.28) 

becomes 

w(t) = exp \-m(t) + tr (ln(l + e- M(t) e- L{t) )^ . (B.29) 

Using Eqs. (B.12), (B.13), (B.16) and (B.17), we get 

w(i) - y(t) exp [tr (ln(l - p(t) - a(t) + 2a(t)p(t)))] . 

(B.30) 

The action contains a term with a time derivative, like 
Tr[i^] in Eq. (2.12), which, starting from Eq. (B.30), 
becomes 



Tr 



2 dD 








= lo tr 


( 


dt 







2a- 1 



dp 

1 — p — a + 2(7 p ) dt 



(B.31) 



The other terms of the action contain the Hamiltonian 
H. For mean-field states, all the information on the sys- 
tem is contained in the one-body density matrix p and we 
can write the expectation value of the Hamiltonian as an 
energy density functional, i.e., 

Tr \b(t)H~\ = E[p(t)}. (B.32) 

Similarly, we have 

Tr \b{t)A{t)H\ = u)(t)E{p'(t)] and (B.33) 
Tr\A(t)D(t)H\ = uj{t)E[a'{t)]. (B.34) 

As a result, the action defined in Eq. (2.12) becomes 
r-ti 



J = u{ti) 



dtw{t) iE[p' (t)] - iE[a' {t)} 



tr 



2<r(f) - 1 



dp(t) \ 



1 - p(t) - a{t) + 2a(t)p(t) dt J 



(B.35) 
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Equations of motion for p, p' and a' 

We now seek for an equation of motion for p(t). This is 
obtained by requiring the stationarity of J when y and a 
vary. J depends on y due to u which is linear in y, i.e., 



duj 
dy 



to 

y ' 



(B.36) 



Requiring the stationarity of the action when y(t) varies 
implies 



= iE[p'(t)} - iE[a'(t)] + tr ( r](t) 



d P (t) \ 



dt I 



(B.37) 



where we have introduced a new matrix 
2a(t) - 1 



1 - p(t) - a(t) + 2a{t)p{t) 
A variation of a implies a variation of r\ according to 



(B.38) 



Srj = -^-5a. 
da 



(B.39) 



Then, instead of considering the variation of a, we equiv- 
alently study the variation of J with rj. The rj matrix can 
be re-written in such a way that it contains an explicit 
dependence on p' and a': 



v = p- 1 ( P ' -pXi-p)- 1 

= (l-p)-\a' -p)p-\ 



(B.40) 
(B.41) 

This tells us how p' and a' vary with a variation of 77: 



5p = pSrj (1 — p) and 
da' = (1 — p) Srj p. 



(B.42) 
(B.43) 



These relations are used to determine the variation of the 
energies E[p'] and £7 [a 7 ]. Note that the latter can be writ- 
ten as a function of the self-consistent HF Hamiltonian h 
using Eq. (A. 9): 



E[p'} = tr (h[p']p r ) and 
E[a'\ =tr(%V). 



(B.44) 
(B.45) 



A variation of r} induces, then, a variation of Eq. (B.37) 
which can be written as 



tr 



[^+i(l-p)h[p'}p-iph[a'}(l-p) 



0. 



(B.46) 

Requesting this equation to hold for any variation br\ gives 
the equation of motion for p: 



. dp 



i— = (l-p) h[p] p - ph[a] (1 - p). 



(B.47) 



To get the equations of motion for p' and a' , it is con- 
venient to express the time derivative term with the L and 
M matrices. From Eq. (B.13), we have 



-M 



I dp I 



dt 



1 — p dt 1 — p 



(B.48) 



Using, Eq. (B.47), we get 
. d 



t e M = h[p'}e- M -e~ M h[a']. 



Similarly, we have 
d 



ij t e~ L = h[a'\e- L - e~ L h[p'}. 



(B.49) 



(B.50) 



We can then write 



d , dp~ M dp~ L 

dt dt dt 



e- M ,h[p'] 



From Eq. (B.26), we have 



1 + e 



-M' 



This leads to 



d de' M ' 



(B.51) 



(B.52) 



(B.53) 



Using Eq. (B.51) we get the equation of motion for p' 
.dp' 



Similarly, we have 



dt 



(B.54) 



(B.55) 



Development in powers of e 

From the definition of uj(€) in Eq. (B.28), we see that the 
expectation of A\ at the final time t\ is equal to uj(ti). We 
now show that, when the action is stationary, uj(t) is in fact 
constant in time. The one-body density matrix associated 
to a Slater determinant has eigenvalues 1 for the occu- 
pied states and for the others. As a result, we see from 
Eq. (C.15) that m(t) is equal to an infinite constant and, 
then, dm/dt = 0. Using this property and Eqs. (B.23), 
(B.26), (B.29), and (B.51), we get 



d d ( 
—uj cx — tr mil 
dt dt \ 



-M 



d_ p -M' 
-A/' 



cx tr 

1 + e- 
cx tr [h[p'],p'] = 0. 



(B.56) 



uj{t) is then a constant we need to evaluate to get the 
fluctuations and correlations from Eq. (2.29). 

lnw can be expressed as a function of p and L from 
Eq. (B.29) and using Eqs. (B.12) and (B.13): 



lnw = tr [in (l - p + pe~ L )] 



(B.57) 
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Deriving In a; according to the parameter £j, we get 



— lnw = —{Qi 



= tr 



p'e L 



j 

der L ~ 



(B.58) 



(B.59) 



where the first identity is obtained from Eq. (2.29) and the 
second from Eqs. (B.57), (B.9) and (B.26). Let us expand 
p and L in powers of the e: 

p(t) = p&> (t) + p (1) (*) + ••• and (B.60) 
L(t) = L^(t) + --- (B.61) 

where the superscript denotes the power in s. The fact 
that = is due to the boundary condition in Eq. (B.2) 
implying 

i 

We now expand p 1 up to first order in e using Eqs. (B.26) 
and (B.13): 



p' = p(°> + - p^L^(l - p(°>) + 0(e 2 ). (B.63) 

;he 

p'w we get 



Separating the zeroth and first order contributions, i.e., 



Expectation value of one-body observables 

At the zeroth order, Eq. (B.59) becomes 



(Qi) 



tr 

tr[p (0) (*i)Q< . 



(B.65) 



where we have used the boundary condition in Eq. (B.62). 
From Eqs. (B.54) and (B.64) we see that p^ obeys the 
TDHF equation. Eq. (B.65) is then exactly the TDHF 
result for the expectation value of one-body observables. 



Fluctuations and correlations 

We now expand Eq. (B.59) up to first order in e, using 
Eqs. (B.64) and (B.62): 

din io 



de.; 



tr 



P (0) (t i ) + (t i ) - p (0) (ti ) £ ei Qj(l - p (0) (*0) 



where we assumed [Qi,<3j] = for simplicity. Identifying 
with Eq. (B.58), we get 

X = E^ tr [QiP (0) (*i)Q^ (i - p (0) (ti) 



tr(pW(*i)Q, 



(B.67) 



The first term in the right hand side gives the usual fluctu- 
ations and correlations computed at the TDHF level [see 
Eq. (2.37)]. The present approach, which is more general 
as it uses a larger variational space for the evolution of 
the observable than in the TDHF formalism, is expected 
to improve the fluctuations and correlations thanks to the 
second term in the right hand side. 

We see in Eq. (B.67) that the time evolution of p^> 
is needed in addition to p(°K Unfortunately, p^> does not 
obey a simple equation and its determination is not trivial. 
It is however possible to find an equation for dj which 
is easier to implement. In fact, although Eq. (B.67) has 
been obtained by evaluating the matrices at time t±, it is 
possible to show that a similar equation holds at any time. 
Let us first introduce the quantity 



F ij (t)= 1 -t T (p^[4 1 \L^]) + 



tr 



tr IpWL 



(1) r(l) 



L^pWLf\l-pW) 
(B.68) 



where the matrices in the right hand side are expressed at 



(B.64) time t and the L\ obey the boundary condition 

4 1) (t 1 )=s i Q i . 
We have also the boundary condition 



(B.69) 



(B.70) 



Note that a similar term than the first term in the right 
hand side of Eq. (B.68) appears in Eq. (B.67) when we 
consider the general case [Qi,Qj] ^ 0. We now show that, 
in fact, Fij does not depend on time. We first re-write F^j 
as a function of the matrices L, a' and p': 



F; 



-tr 



Pi 



(B.71) 



where we have used Eq. (B.64) and its equivalent for a'. 
The matrices p'' 1 - 1 and a'^ are obtained by taking the first 
order in e in the equations (B.54) and (B.55), respectively. 
We obtain the time-dependent RPA equations 



l : 

dt 
da'W 

! : 

dt 



/i[p (0) ]y (1) l + \tr 2 (vp'^),M (B.72) 



h[p 



(0)1 ^(i) 



tT 2 (va'W),pW] ,(B.73) 



where v is the antisymetrised two-body interaction defined 
in Eq. (2.21). We also get a similar equation for L: 



1 + E SnQn -Qi + E £ m QiQr. 



(B.66) 



l : 

dt 



h\p®],LW 



tr 2 



v[LW,pW] 



(B.74) 
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The time derivative of F^j obeys 
1 



dt 13 



-tr 



dL[ 1} , (1) 



dt 



H dt 



-d) 



da 1 



(i) 



dt 



(B.75) 



Replacing with Eqs. (B.72), (B.73) and (B.74), we get 

d 
dt 



-F ij {t) = Q. 



(B.76) 



Fij is then constant with time. 

It is convenient to write it with the matrices at the time 
to instead of t\. Indeed, in this case, we have p (1) (t ) = 
due to the boundary condition at initial time. This sim- 
plifies the expression for the fluctuations and correlations: 



L i 



lim 



1 



>0 2£i£j 



tr 



[Lf\t ),p^{to)][Lf\to),p { - Q \to)]} (B.77) 



Let us introduce a new matrix r)i(t,Si) defined as 



Rje ilW(*) p (0) (t)e -iLW(t) | 



(B.78) 
(B.79) 



where the last equation holds for small e. We can show 
that rii(t, Ei) follows also a TDHF equation. Note that the 
fact that is a Slater determinant implies that rji is also 
a Slater determinant according to Eq. (B.79). 

From Eqs. (B.77) and (B.78) we get our final result 



lim 



tr 



with the boundary condition at the final time t\ 



r 1i {t 1 ,e l )=e^p^\t 1 )e-^ 



(B.80) 



(B.81) 



C From Fock to single-particle space with 
exponential of one-body operators 

Several relations are derived which are used in the deriva- 
tion of the Balian-Veneroni prescription for fluctuation 
and correlation of one-body operators in appendix B. They 
relate operators or their associated matrices in Fock space 
to matrices expressed in the single-particle space. 



Link between the density matrix of an independent- 
particle state and the one-body density matrix 

Our goal is to show the relation 



1 + e 



A I 



(C.l) 



where p is the one-body density matrix with elements 

Pap = (a^a a ) = Tr {jja^do^j , (C.2) 

and 

D = e - m - M (C.3) 

is the density matrix associated to an independent-particle 
state described by a Slater determinant. The later obeys 

Tib = f (C.4) 

and M is a one-body operator of the form 

M = Y,M a paid p . (C.5) 

a/3 



We first show the following property: 



(C.6) 



Define F(x) = e Be . Its Taylor development reads 



x n d n F 
F( X ) = £ -/„ with f n = — 

n—0 



Tl = 



B + x A B 



m=0 

x 2 
+ 2! 



[A, [A, b] 



=>■ fa = B, fn+i = [A, /„]. Using this relation with x = 1 , 
we get 



oo 

D- l a^D = V -/„ with 



Tl = 



fn+l = M, f n 



where the f n have an implicit label a. Using Eq. (C.5) 
and the commutation rules for creators and annihilators 
of fcrmions, we find 

;„=^(M")^4. 

(8 



We finally get 



leading to Eq. (C.6). 

Using a) a ap = 5 a p — d/sat, and Eqs. (C.2) and (C.6) 
leads directly to 



a = 8f) a - Tr (dlbd/3 J = 6f) a - ^ {e M ) Pfi T (C.7) 



Pl3c 



We get p = 1 — pe and finally Eq. (C.l) 
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Trace of the exponential of a one-body operator 

Our goal is to show the relation 



Tre 



-A 



det(l + e- A ) 



where A is a one-body operator. 
We first show that 

TrB = [ dz*dz e~ z 



(z\B\-z) 



(C.9) 



where B contains only even powers of and & (this is the 
case, for instance, when B is the exponential of a one-body 
operator), and z is an element of the Grassmann algebra 
(see appendix D). 

With {|0} a basis of the Fock space, the trace of B 
reads TrB = Inserting the closure relation 

from Eq. (D.ll), we get 



TrB = J dz*dze- z * z J2^\B\z)(z\0- 



(CIO) 



From Eq. (D.12), we have 



{|£)} can be chosen to be an ensemble of Slater determi- 
nants with different particle numbers N^: 

|0=(IK 

Using Eq. (D.3) we get 

m= ilk 

Indeed, from Eq. (D.4), we have 



N. 



I-) = 1. 



Then, using Eq. (D.2), we get 

(t\z) (z\Z} = (-l) N t(z\Z)(Z\z) = {z\t}(t\-z). 

Using the closure relation J2 10(0 = 1> Eq. (CIO) finally 
becomes 

TtB = [ dz*dz e~ z * z B(V Z , z) (z\ - z) 
dz*dz e~ z * z (z\ 

i.e., Eq. (C.9). 

The second step is to show the relation 



e- A \z) = \e- A z) 



(C.ll) 



Similarly to Eq. (C.6), we have 



^- A = e- A J2{e A ) nfi ap 



(C-8) an d then, using Eq. (D.3), we get 



—A 



E 



'a/3 



(e~ A z) e~ A \z) = a a \e~ A z) 



Finally, we write Eq. (C.9) with B = e A and we use 
Eqs. (C.ll) and (D.10) to get 

Tre- 4 = J dz*dz e~ z * z {z\ - e~ A z) 

= J dz*dz e - z '( 1+e ~ A >. 

Using Eq. (D.14) we obtain our final result in Eq. (C.8). 

Trace of the density matrix 

Our goal is to show the relation 



TrD = e -™+tr[ln(l+e- M )] 



(C.12) 



where Tr denotes the trace in the Fock space while tr de- 
notes the trace of a single-particle matrix. D is the density 
matrix of an independent particle system, as defined by 
Eqs. (C.3) and (C.5). However we do not assume a nor- 
malised state, i.e., TrD is not necessarily equal to 1. 
From Eq. (C.8), we have 

Trr) = e- m det(l + e - M )- (C.13) 

We now show the property 

trln(l + A) = lndet(l + A) (C.14) 

where A is a matrix which can be diagonalised. In its 
diagonal form, Ay = a,<%, we have trA' 1 = J^i a ?- Then 
we get 

trln(l + A) = tr(A - - A 2 + - 



i 

= + a 

i 

= lnIJ(l + Oi) 

i 

= lndet(l + A). 



A A - 

i 

1 , 



) 



(C.15) 

(C.16) 

- M and 



Taking the exponential of Eq. (C.14) with A 
replacing in Eq. (C.13) gives Eq. (C.12). 

In the particular case where the system is described by 
a Slater determinant, the state is normalised (TrD = 1) 
and one gets 



m = tr ln(l 



) 



(C.17) 
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D Grassmann algebra 



A brief introduction to Grassmann algebra can be found 
in p. 21-27 of Rcf. [2]. Here, we recall mostly the main 
relations which are used in appendix C. Similarly to the 
imaginary number i which has been introduced to satisfy 
i 2 = — 1, the Grassmann algebra has been developed such 
that its elements vanish when squared. Elements {z a } of 
the Grassmann algebra and their complex conjugated el- 
ements obey 



2 — * 2 — n 

z a — z a ~ U 



(D.l) 



Z a Zp + Zf3Z a = Z a Z*p + ZpZ a = 0. (D.2) 

They are eigenvalues of annihilation and creation opera- 
tors such that 

a a \z) = z a \z) and {z\a) a = {z\z* a (D.3) 

where \z) and (z\ are coherent states of the form 

\ z ) = e £«,*.a£|_) an d ( z \ = (-| e E- , (D.4) 

and |— ) is the particle vacuum. 

Elements of the Grassmann algebra are polynomial of 
second degree at most, as, for example, 

P(z a> z* a ) = P + Piz a + P 2 z* + P 12 z a z* a . (D.5) 

Derivation is defined as 

t Za P = P 1 +P 12 z* a 

~3 z , a p = p 2 - p 12Za 

(indeed ~c) z * z a z* a = -(d z *z*)z a = -z a ) 



d Za P = Pi - Pl2Z* a 



1— 



d Z ,P = P 2 + Pl2Z a . 



We can show the following properties 



(D.6) 



il\z) = d z Jz) and (z\a a = d z *(z\. (D.7) 



Integration is defined as 



J dz a = J dz* a = 



and 



dz q, -2q, — 1 . 



J dz a z a - 

The overlap of two coherent states reads 



(z\z>) = e 



(D. 



(D.9) 



(D.10) 



Indeed, using Eqs. (D.l) and (D.4), and the Wick theorem, 
we have 

n \ a / rn \ /3 J 

= 1 + ^2 z * a (-\a a a\\-)z'p 



a/3 

-y 

2! 2 ^ 



Za^{-\a a a p a\al\-)zi f z' s + ' ' ' 



1_ <~ Z a Z a ' 2?2 2-^^ ZaZ P Z P Za ~ Z a Z /3 Z a Z p) + ' ' ' 



a/3 



1 + 2^ z * z <* + z a z * 1^ z p z p i~ " ' 

a a ft 



Eq. (D.10) can be used to define the metric entering the 
following closure relation 



dz*dz e" z z \z) (z\ = 1, 



(D.ll) 



where we introduced the notation J dz*dz = JYl a dz^dz a . 

An operator A{a\ a) = A(- ■ ■ , , ■ ■ • , dp, • ■ ■ ) may be 
represented by the differential operator A defined as 



(z\A{a\a)\z')=A(d*,z'){z\z') 



(z\A(a?,a)\z')=A(z*, d z .)(z\z') 



(D.12) 



(D.13) 



where we used the notation d z = {• • • , d Za , • • • } and 

9 X = {--- ,o Za ,---}- 

We now show the relation 



dz*dz e~ z * Mz = det M (D.14) 



Starting from the Taylor development 



1 

e~ z * Mz = / dz*dz V — (-z*Mz) n , 
/ n! 
17 n=0 



and using Eq. (D.8), we see that the only non- vanishing 
term of the Taylor development is for n = N where N is 
the dimension of the Fock space: 



dz* dz e 



(-1) 



AT 



AH 



dz*dz (z*Mz) 



N 



with 



(Z*MZ) N = 4 M «lA^ft '" E Z a« M "«Av^ 

a 1 /3 1 
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From Eq. (D.l), we see that all terms with = OLj^i or 
1% = fij^i vanish. It leads to two sums over all possible 
permutations of the a and of the /3 indices: 

(Z*MZ) N = <(l) M a(l)/3(l) Z /3(l) • • • 

a,/3£Perm jy 

Z a(N) M a(N)0(N)Zf3(N)- 

Noting the following relations: 

J dz*dz = (-l)N(N-l)/2 

J dz* ■ ■ ■ dz* N dz\ ■ ■ ■ dzN, 

Z a(l) Z l3(l)- ■ ■ Z a(N) Z /3(N) = (-1) 

2/3(1) ■ ■ ■ z /3(N) z a(l) ■ ■ ■ z u{N)i 

z m ■ ■ ■ z p { N) = (-l) N{N - 1),2 smn{P)z N ...a*, 
and the similar equation for the z*, we get 

J dz*dz z a(1) *zp (1) ■ ■ ■ z a{N) *z l3{N) = (-l) Jv sign(a)sign(/3) 

where sign(a) is the sign of the permutation a, i.e., +1 for 
an even number of permutations and -1 for an odd one. 
This leads to Eq. (D.14) 

J ' dz*dze- z ' Mz = i Y, sign(a)sign(/3) 

M a (l)/3(1) ■ ■ ■ M a (N)fi(N) 
= Y si S n ( a ) M a(l)l ■ ■ ■ M a ( N ) N 

= det M. 
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